hera_sim¶

Basic simulation package for HERA-like redundant interferometric arrays.

Features¶

Systematic Models: Many models of instrumental systematics in various forms, eg. thermal noise, RFI, bandpass gains, cross-talk, cable reflections and foregrounds.
HERA-tuned: All models have defaults tuned to HERA, with various default “sets” available (eg.H1C, H2C)
Interoperability: Interoperability with pyuvdata datasets and pyuvsim configurations.
Ease-of-use: High-level interface for adding multiple systematics to existing visibilities in a self-consistent way.
Visibility Simulation: A high-level interface for visbility simulation that is compatible with the configuration definition from pyuvsim but is able to call multiple simulator implementations.
Convenience: Methods for adjusting simulated data to match the times/baselines of a reference dataset.

Documentation¶

At ReadTheDocs. In particular, for a tutorial and overview of available features, check out the tour.

Installation¶

Conda users¶

If you are using conda, the following command will install all dependencies which it can handle natively:

$ conda install -c conda-forge numpy scipy pyuvdata attrs h5py healpy pyyaml

If you are creating a new development environment, consider using the included environment file:

$ conda env create -f ci/tests.yaml

This will create a fresh environment with all required dependencies, as well as those required for testing. Then follow the pip-only instructions below to install hera_sim itself.

Pip-only install¶

Simply use pip install -e . or run pip install git+git://github.com/HERA-Team/hera_sim.

Developer install¶

For a development install (tests and documentation), run pip install -e .[dev].

Other optional extras can be installed as well. To use baseline-dependent averaging functionality, install the extra [bda]. For the ability to simulate redundant gains, install [cal]. To enable GPU functionality on some of the methods (especially visibility simulators), install [gpu].

As the repository is becoming quite large, you may also wish to perform a shallow clone to retrieve only the recent commits and history. This makes the clone faster and avoid bottleneck in CI pipelines.

Provide an argument --depth 1 to the git clone command to copy only the latest revision of the repository.

git clone -–depth [depth] git@github.com:HERA-Team/hera_sim.git

Versioning¶

We use semantic versioning (major.minor.patch) for the hera_sim package (see SemVer documentation). To briefly summarize, new major versions include API-breaking changes, new minor versions add new features in a backwards-compatible way, and new patch versions implement backwards-compatible bug fixes.

Contents¶

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.

[1]:

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

[2]:

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)

[3]:

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¶

[4]:

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

[5]:

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¶

[6]:

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

[7]:

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¶

[8]:

vis_fg = vis_fg_diffuse + vis_fg_pntsrc

[9]:

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¶

[10]:

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)

[11]:

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()

[12]:

vis_fg_nos = vis_fg + nos_jy

[13]:

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¶

[14]:

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

[15]:

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()

[16]:

vis_fg_nos_rfi = vis_fg_nos + rfi_all

[17]:

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¶

[18]:

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)

[19]:

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¶

[20]:

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).

This notebook provides a brief tutorial on basic use of the Simulator class, followed by a longer, more in-depth tutorial that shows some advanced features of the class.

Setup¶

[1]:

import tempfile
import time
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
from astropy import units

import hera_sim
from hera_sim import Simulator, DATA_PATH, utils
from uvtools.plot import labeled_waterfall

We’ll inspect the visibilities as we go along by plotting the amplitudes and phases on different axes in the same figure. Here’s the function we’ll be using:

[2]:

def waterfall(sim, antpairpol=(0,1,"xx"), figsize=(6,3.5), dpi=200, title=None):
    """Convenient plotting function to show amp/phase."""
    fig, (ax1, ax2) = plt.subplots(
        nrows=2,
        ncols=1,
        figsize=figsize,
        dpi=dpi,
    )
    fig, ax1 = labeled_waterfall(
        sim.data,
        antpairpol=antpairpol,
        mode="log",
        ax=ax1,
        set_title=title,
    )
    ax1.set_xlabel(None)
    ax1.set_xticklabels(['' for tick in ax1.get_xticks()])
    fig, ax2 = labeled_waterfall(
        sim.data,
        antpairpol=antpairpol,
        mode="phs",
        ax=ax2,
        set_title=False,
    )
    return fig

Introduction¶

The Simulator class provides an easy-to-use interface for all of the features provided by hera_sim, though it is not to be confused with the VisibilitySimulator class, which is intended to be a universal interface for high-accuracy visibility simulators commonly used within the community. The Simulator class features:

A universal add method for simulating any effect included in the hera_sim API.
- More advanced users can create their own custom effects that the Simulator can use, as long as these custom effects abide by a certain syntax.
A get method for retrieving any previously simulated effect.
- Simulated effects are not cached, but rather the parameters characterizing the effect are recorded and the desired effect is re-simulated.
Four modes of specifying how to configure the random state when simulating effects:
- “redundant” ensures that the random state is set the same for each baseline within a redundant group;
- “once” uses the same random state when simulating data for every baseline in the array;
- “initial” only sets the random state initially, prior to simulating data for the array;
- Similar to “once”, the user can provide an integer to seed the random state, which is then used to set the random state for every baseline in the array.
Convenient methods for writing the data to disk.
- The write method writes the entire UVData object’s contents to disk in the desired format.
- The chunk_sim_and_save method allows for writing the data to disk in many files with a set number of integrations per file.

In order to enable a simple, single interface for simulating an effect, hera_sim employs a certain heirarchy to how it thinks about simulated effects. Every simulation effect provided in hera_sim derives from a single abstract base class, the SimulationComponent. A “component” is a general term for a category of simulated effect, while a “model” is an implementation of a particular effect. Components track all of the models that are particular implementations of the components, and this enables hera_sim to track every model that is created, so long as it is a subclass of a simulation component. To see a nicely formatted list of all of the components and their models, you can print the output of the components.list_all_components function, as demonstrated below.

[3]:

print(hera_sim.components.list_all_components())

array:
  lineararray
  hexarray
foreground:
  diffuseforeground | diffuse_foreground
  pointsourceforeground | pntsrc_foreground
noise:
  thermalnoise | thermal_noise
rfi:
  stations | rfi_stations
  impulse | rfi_impulse
  scatter | rfi_scatter
  dtv | rfi_dtv
gain:
  bandpass | gains | bandpass_gain
  reflections | reflection_gains | sigchain_reflections
crosstalk:
  crosscouplingcrosstalk | cross_coupling_xtalk
  crosscouplingspectrum | cross_coupling_spectrum | xtalk_spectrum
  whitenoisecrosstalk | whitenoise_xtalk | white_noise_xtalk
eor:
  noiselikeeor | noiselike_eor

The format for the output is structured as follows:

component_1:
  model_1_name | model_1_alias_1 | ... | model_1_alias_N
  ...
  model_N_name | model_N_alias_1 | ...
component_2:
  ...

So long as there are not name conflicts, hera_sim is able to uniquely identify which model corresponds to which name/alias, and so there is some flexibility in telling the Simulator how to find and simulate an effect for a particular model.

Basic Use¶

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 using an integration time of 10.7 seconds, and a 7-element hexagonal array.

[4]:

# Define the array layout.
array_layout = hera_sim.antpos.hex_array(
    2, split_core=False, outriggers=0
)

# Define the timing parameters.
start_time = 2458115.9  # JD
integration_time = 10.7  # seconds
Ntimes = int(30 * units.min.to("s") / integration_time)

# Define the frequency parameters.
Nfreqs = 100
bandwidth = 1e8  # Hz
start_freq = 1e8  # Hz

sim_params = dict(
    Nfreqs=Nfreqs,
    start_freq=start_freq,
    bandwidth=bandwidth,
    Ntimes=Ntimes,
    start_time=start_time,
    integration_time=integration_time,
    array_layout=array_layout,
)

# Create an instance of the Simulator class.
sim = Simulator(**sim_params)

Overview¶

The Simulator class adds some attributes for conveniently accessing metadata:

[5]:

# Observed frequencies in GHz
sim.freqs[::10]

[5]:

array([0.1 , 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19])

[6]:

# Observed Local Sidereal Times (LSTs) in radians
sim.lsts[::10]

[6]:

array([4.58108965, 4.58889222, 4.59669478, 4.60449734, 4.61229991,
       4.62010247, 4.62790503, 4.6357076 , 4.64351016, 4.65131273,
       4.65911529, 4.66691785, 4.67472042, 4.68252298, 4.69032555,
       4.69812811, 4.70593067])

[7]:

# Array layout in local East-North-Up (ENU) coordinates
sim.antpos

[7]:

{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])}

[8]:

# Polarization array
sim.pols

[8]:

['xx']

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

[9]:

fig = sim.plot_array()

_images/tutorials_hera_sim_simulator_21_0.png

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

[10]:

type(sim.data)

[10]:

pyuvdata.uvdata.uvdata.UVData

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 an appropriately defined callable class (see the section on using custom classes for exact details), while the variable keyword arguments parametrize the model. The add method simulates the effect and applies it to the entire array. Let’s walk through an example.

We’ll start by simulating diffuse foreground emission. For simplicity, we’ll be using the H1C season default settings so that the number of extra parameters we need to specify is minimal.

[11]:

# Use H1C season defaults.
hera_sim.defaults.set("h1c")

[12]:

# Start off by simulating some foreground emission.
sim.add("diffuse_foreground")

You have not specified how to seed the random state. This effect might not be exactly recoverable.

[13]:

# Let's check out the data for the (0,1) baseline.
fig = waterfall(sim, title="Simulated Diffuse Foregrounds")
fig.tight_layout()

FixedFormatter should only be used together with FixedLocator

_images/tutorials_hera_sim_simulator_29_1.png

That was simple enough; however, basic use like this does not ensure that the simulation is as realistic as it can be with the tools provided by hera_sim. For example, data that should be redundant is not redundant by default:

[14]:

# (0,1) and (5,6) are redundant baselines, but...
np.all(sim.get_data(0,1,"xx") == sim.get_data(5,6,"xx"))

[14]:

False

[15]:

# As an extra comparison, let's plot them
fig1 = waterfall(sim, antpairpol=(0,1,"xx"), title="Baseline (0,1)")
fig2 = waterfall(sim, antpairpol=(5,6,"xx"), title="Baseline (5,6)")

FixedFormatter should only be used together with FixedLocator

_images/tutorials_hera_sim_simulator_32_1.png

_images/tutorials_hera_sim_simulator_32_2.png

We can ensure that a simulated effect that should be redundant is redundant by specifying that we use a “redundant” seed:

[16]:

sim.refresh()  # Clear the contents and zero the data
sim.add(
    "diffuse_foreground",
    seed="redundant",  # Use the same random state for each redundant group
)

[17]:

np.all(sim.get_data(0,1,"xx") == sim.get_data(5,6,"xx"))

[17]:

True

[18]:

# Let's plot them again in case the direct comparison wasn't enough.
fig1 = waterfall(sim, antpairpol=(0,1,"xx"), title="Baseline (0,1)")
fig2 = waterfall(sim, antpairpol=(5,6,"xx"), title="Baseline (5,6)")

FixedFormatter should only be used together with FixedLocator

_images/tutorials_hera_sim_simulator_36_1.png

_images/tutorials_hera_sim_simulator_36_2.png

It is worth noting that the get_data method used here is the pyuvdata.UVData method that has been exposed to the Simulator, not to be confused with the get method of the Simulator, but more on that later.

[19]:

# Go one step further with this example.
all_redundant = True
for red_grp in sim.red_grps:
    if len(red_grp) == 1:
        continue  # Nothing to see here
    base_vis = sim.get_data(red_grp[0])
    # Really check that data is redundant.
    for baseline in red_grp[1::]:
        all_redundant &= np.allclose(base_vis, sim.get_data(baseline))
all_redundant

[19]:

True

Retrieving Simulated Effects¶

The parameter values used for simulating an effect are stored so that the effect may be recovered at a later time by re-simulating the effect. For effects that have a random element (like bandpass gains or thermal noise), you’ll want to make sure that the seed parameter is specified to ensure that the random state can be configured appropriately.

[20]:

sim.refresh()
vis = sim.add(
    "noiselike_eor",
    eor_amp=1e-3,
    seed="redundant",
    ret_vis=True,  # Return a copy of the visibility for comparison later
)
sim.add(
    "gains",  # Apply a bandpass to the EoR data
    seed="once",  # Use the same seed each pass through the loop
)

As an extra note regarding the "once" setting used to set the random state in the above cell, this setting ensures that the same random state is used each iteration of the loop over the array. This isn’t strictly necessary to use for multiplicative effects, since these are computed once before iterating over the array, but should be used for something like point-source foregrounds (which randomly populates the sky with sources each time the model is called).

[21]:

# The data has been modified by the gains, so it won't match the original.
np.allclose(vis, sim.data.data_array)

[21]:

False

[22]:

# We can recover the simulated effect if we wanted to, though.
np.allclose(vis, sim.get("noiselike_eor"))

[22]:

True

It is worth highlighting that the get method does not retrieve cached data (since simulated effects are not cached), but rather re-simulates the effect using the same parameter values and random states (since these are cached).

[23]:

# Make explicitly clear that the simulated effects are not cached
vis is sim.get("noiselike_eor")

[23]:

False

[24]:

# Show that the same effect really is recovered.
gains = sim.get("gains")
new_vis = vis.copy()
for ant1, ant2, pol, blt_inds, pol_ind in sim._iterate_antpair_pols():
    # Modify the old visibilities by the complex gains manually.
    total_gain = gains[(ant1, pol[0])] * np.conj(gains[(ant2, pol[1])])
    new_vis[blt_inds,0,:,pol_ind] *= total_gain
np.allclose(new_vis, sim.data.data_array)  # Now check that they match.

[24]:

True

Each time a component is added to the simulation, it is logged in the history:

[25]:

print(sim.data.history)

hera_sim v1.0.1.dev14+ga590958: Added noiselikeeor using parameters:
eor_amp = 0.001
fringe_filter_type = tophat
hera_sim v1.0.1.dev14+ga590958: Added bandpass using parameters:
bp_poly = <hera_sim.interpolators.Bandpass object at 0x7f124d22d550>

Saving A Simulation¶

Finally, we’ll often want to write simulation data to disk. The simplest way to do this is with the write method:

[26]:

tempdir = Path(tempfile.mkdtemp())
filename = tempdir / "simple_example.uvh5"
sim.write(filename, save_format="uvh5")

[27]:

filename in tempdir.iterdir()

[27]:

True

[28]:

# Check that the data is the same.
sim2 = Simulator(data=filename)
is_equiv = True
# Just do a basic check.
for attr in ("data_array", "freq_array", "time_array", "antenna_positions"):
    is_equiv &= np.all(getattr(sim.data, attr) == getattr(sim2.data, attr))
is_equiv

Telescope hera_sim is not in known_telescopes.

[28]:

True

This concludes the section on basic use of the Simulator.

Advanced Use¶

The preceding examples should provide enough information to get you started with the Simulator. The rest of this notebook shows some of the more advanced features offered by the Simulator.

Singling Out Baselines/Antennas/Polarizations¶

It is possible to simulate an effect for only a subset of antennas, baselines, and/or polarizations by making use of the vis_filter parameter when using the add method. Below are some examples.

[29]:

sim.refresh()
vis_filter = [(0,1), (3,6)]  # Only do two baselines
sim.add("diffuse_foreground", vis_filter=vis_filter)
for baseline in sim.data.get_antpairs():
    if np.any(sim.get_data(baseline)):
        print(f"Baseline {baseline} has had an effect applied.")

Baseline (1, 0) has had an effect applied.
Baseline (6, 3) has had an effect applied.

You have not specified how to seed the random state. This effect might not be exactly recoverable.

[30]:

vis_filter = [1,]  # Apply gains to only one antenna
vis_A = sim.get_data(0, 1)
vis_B = sim.get_data(3, 6)
sim.add("gains", vis_filter=vis_filter)
np.allclose(vis_A, sim.get_data(0, 1)), np.allclose(vis_B, sim.get_data(3, 6))

[30]:

(False, True)

[31]:

# Now something a little more complicated.
hera_sim.defaults.deactivate()  # Use the same initial params, but 2 pol
sim = Simulator(polarization_array=["xx", "yy"], **sim_params)
hera_sim.defaults.activate()

# Actually calculate the effects for the example.
vis_filter = [(0,1,"yy"), (0,2,"yy"), (1,2,"yy")]
sim.add("noiselike_eor", vis_filter=vis_filter)
for antpairpol, vis in sim.data.antpairpol_iter():
    if np.any(vis):
        print(f"Antpairpol {antpairpol} has had an effect applied.")
# Apply gains only to antenna 0
vis_A = sim.get_data(0, 1, "yy")
vis_B = sim.get_data(0, 2, "yy")
vis_C = sim.get_data(1, 2, "yy")
sim.add("gains", vis_filter=[0,])
(
    np.allclose(vis_A, sim.get_data(0, 1, "yy")),
    np.allclose(vis_B, sim.get_data(0, 2, "yy")),
    np.allclose(vis_C, sim.get_data(1, 2, "yy")),
)

You have not specified how to seed the random state. This effect might not be exactly recoverable.

Antpairpol (1, 0, 'yy') has had an effect applied.
Antpairpol (2, 0, 'yy') has had an effect applied.
Antpairpol (2, 1, 'yy') has had an effect applied.

[31]:

(False, False, True)

For some insight into what’s going on under the hood, this feature is implemented by recursively checking each entry in vis_filter and seeing if the current baseline + polarization (in the loop that simulates the effect for every baseline + polarization) matches any of the keys provided in vis_filter. If any of the keys are a match, then the filter says to simulate the effect—this is important to note because it can lead to some unexpected consequences. For example:

[32]:

seed = 12345  # Ensure that the random components are identical.
sim.refresh()
vis_filter = [(0,1), "yy"]
vis_A = sim.add(
    "diffuse_foreground",
    vis_filter=vis_filter,
    ret_vis=True,
    seed=seed,
)
sim.refresh()
vis_filter = [(0,1,"yy"),]
vis_B = sim.add(
    "diffuse_foreground",
    vis_filter=vis_filter,
    ret_vis=True,
    seed=seed,
)
np.allclose(vis_A, vis_B)

[32]:

False

[33]:

# In case A, data was simulated for both polarizations for baseline (0,1)
blt_inds = sim.data.antpair2ind(0, 1, ordered=False)
np.all(vis_A[blt_inds])

[33]:

True

[34]:

# As well as every baseline for polarization "yy"
np.all(vis_A[...,1])

[34]:

True

[35]:

# Only baseline (0,1) had an effect simulated for polarization "xx"
np.all(vis_A[...,0])

[35]:

False

[36]:

# Whereas for case B, only baseline (0, 1) with polarization "yy"
# had the simulated effect applied.
(
    np.all(vis_B[blt_inds,...,1]),  # Data for (0, 1, "yy")
    np.all(vis_B[blt_inds]),  # Data for both pols of (0, 1)
    np.all(vis_B[...,1]),  # Data for all baselines with pol "yy"
)

[36]:

(True, False, False)

Here’s an example of how the vis_filter parameter can be used to simulate a single antenna that is especially noisy.

[37]:

hera_sim.defaults.deactivate()  # Pause the H1C defaults for a moment...
sim = Simulator(**sim_params)  # Only need one polarization for this example.
hera_sim.defaults.activate()  # Turn the defaults back on for simulating effects.
# Start by adding some foregrounds.
sim.add("diffuse_foreground", seed="redundant")
# Then add noise with a receiver temperature of 100 K to all baselines
# except those that contain antenna 0.
v = sim.get_data(0, 1)
vis_filter = [antpair for antpair in sim.get_antpairs() if 0 not in antpair]
sim.add("thermal_noise", Trx=100, seed="initial", vis_filter=vis_filter)
# Now make antenna 0 dramatically noisy.
sim.add(
    "thermal_noise",
    Trx=5e4,
    seed="initial",
    vis_filter=[0,],
    component_name="noisy_ant",
)

[38]:

# Recall that baselines (0,1) and (2,3) are redundant,
# so the only difference here is the noise.
jansky_to_kelvin = utils.jansky_to_kelvin(
    sim.freqs,
    hera_sim.defaults('omega_p'),
)
Tsky_A = sim.get_data(0, 1) * jansky_to_kelvin
Tsky_B = sim.get_data(2, 3) * jansky_to_kelvin
np.std(np.abs(Tsky_A)), np.std(np.abs(Tsky_B))

[38]:

(50.151933214628464, 47.48252263186864)

Notice that the final call to the add method specified a value for the parameter component_name. This is especially useful when simulating an effect multiple times using the same class, as we have done here. As currently implemented, the Simulator would normally overwrite the parameters from the first application of noise with the parameters from the second application of noise; however, by specifying a name for the second application of noise, we can recover the two effects independently:

[39]:

vis_A = sim.get("thermal_noise")
vis_B = sim.get("noisy_ant")
blt_inds = sim.data.antpair2ind(0, 1)
(
    np.any(vis_A[blt_inds]),  # The first application give (0, 1) noise,
    np.all(vis_B[blt_inds]),  # but the second application did.
)

[39]:

(False, True)

With that, we conclude this section of the tutorial. This section should have provided you with sufficient examples to add effects to your own simulation in rather complicated ways—we recommend experimenting a bit on your own!

Naming Components¶

This was briefly touched on in the previous section, but it is interesting enough to get its own section. By using the component_name parameter, it is possible to use the same model multiple times for simulating an effect and have the Simulator be able to recover every different effect simulated using that model. Here are some examples:

[40]:

hera_sim.defaults.set("debug")
sim = Simulator(
    Ntimes=6*24,  # Do a full day
    integration_time=10*60,  # In 10 minute increments
    Nfreqs=100,  # Just to make the plots look nicer
)
sim.add(
    "pntsrc_foreground",
    nsrcs=5000,  # A lot of sources
    Smin=0.01,
    Smax=0.5,  # That are relatively faint
    component_name="faint sources",
    seed="once",
)
sim.add(
    "pntsrc_foreground",
    nsrcs=10,  # Just a few sources
    Smin=0.5,
    Smax=10,  # That are bright
    component_name="bright_sources",
    seed="once",
)

[41]:

# Let's look at the data. First, we need to recover it.
faint_srcs = sim.get("faint sources", key=(0,1,"xx"))
bright_srcs = sim.get("bright_sources", key=(0,1,"xx"))

# Now let's get ready to plot.
fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(6,3.5), dpi=200)

# Use a common colorscale
vmin = min(np.abs(faint_srcs).min(), np.abs(bright_srcs).min())
vmax = max(np.abs(faint_srcs).max(), np.abs(bright_srcs).max())

# Actually make the plots
fig, ax1 = labeled_waterfall(
    faint_srcs,
    freqs=sim.freqs * units.GHz.to("Hz"),
    lsts=sim.lsts,
    set_title="Faint Sources",
    ax=ax1,
    mode="abs",
    vmin=vmin,
    vmax=vmax,
)
ax1.xaxis.set_visible(False)
fig, ax2 = labeled_waterfall(
    bright_srcs,
    freqs=sim.freqs * units.GHz.to("Hz"),
    lsts=sim.lsts,
    set_title="Bright Sources",
    mode="abs",
    ax=ax2,
    vmin=vmin,
    vmax=vmax,
)
fig.tight_layout()

_images/tutorials_hera_sim_simulator_78_0.png

So, by setting a value for the component_name parameter when calling the add method, it is possible to use a single model for simulating multiple different effects and later on recover the individual effects.

Pre-Computing Fringe-Rate/Delay Filters¶

Some of the effects that are simulated utilize a fringe-rate filter and/or a delay filter in order to add to the realism of the simulation while keeping the computational overhead relatively low. Typically these filters are calculated on the fly; however, this can start to become an expensive part of the calculation for very large arrays. In order to address this, we implemented the ability to pre-compute these filters and use the cached filters instead. Here’s an example.

[42]:

# Make a big array with a high degree of redundancy.
big_array = hera_sim.antpos.hex_array(6, split_core=False, outriggers=0)
hera_sim.defaults.set("debug")  # Use short time and frequency axes.
sim = Simulator(array_layout=big_array)
delay_filter_kwargs = {"standoff": 30, "delay_filter_type": "gauss"}
fringe_filter_kwargs = {"fringe_filter_type": "gauss", "fr_width": 1e-4}
seed = 12345

[43]:

# First, let's simulate an effect without pre-computing the filters.
vis_A = sim.add(
    "diffuse_foreground",
    delay_filter_kwargs=delay_filter_kwargs,
    fringe_filter_kwargs=fringe_filter_kwargs,
    seed=seed,
    ret_vis=True,
)

[44]:

# Now do it using the pre-computed filters.
sim.calculate_filters(
    delay_filter_kwargs=delay_filter_kwargs,
    fringe_filter_kwargs=fringe_filter_kwargs,
)
vis_B = sim.add("diffuse_foreground", seed=seed, ret_vis=True)

[45]:

# Show that we get the same result.
np.allclose(vis_A, vis_B)

[45]:

True

The above example shows how to use the filter caching mechanism and that it produces results that are identical to those when the filters are calculated on the fly, provided the filters are characterized the same way. Note, however, that this feature is still experimental and extensive benchmarking needs to be done to determine when it is beneficial to use the cached filters.

Using Custom Models¶

In addition to using the models provided by hera_sim, it is possible to make your own models and add an effect to the simulation using the custom model and the Simulator. If you would like to do this, then you need to follow some simple rules for writing the custom model:

A component base class must be defined using the @component decorator.
Models must inherit from the component base class.
The model call signature must have freqs as a parameter. Additive components (like visibilities) must also take lsts as a parameter.
- Additive components must return a np.ndarray with shape (lsts.size, freqs.size).
- Multiplicative components must have a class attribute is_multiplicative that is set to True, and they must return a dictionary mapping antenna numbers to np.ndarrays with shape (freqs.size,) or (lsts.size, freqs.size).
Frequencies will always be passed in units of GHz, and LSTs will always be passed in units of radians.

An example of how to do this is provided in the following section.

Registering Classes¶

[46]:

# Minimal example for using custom classes.
# First, make the base class that the custom model inherits from.
from hera_sim import component

@component
class Example:
    # The base class doesn't *need* a docstring, but writing one is recommended
    pass

class TestAdd(Example):
    def __init__(self):
        pass

    def __call__(self, lsts, freqs):
        return np.ones((len(lsts), len(freqs)), dtype=complex)

sim.refresh()
sim.add(TestAdd)  # You can add the effect by directly referencing the class
np.all(sim.data.data_array == 1)

You have not specified how to seed the random state. This effect might not be exactly recoverable.

[46]:

True

[47]:

# Now for a multiplicative effect
class TestMult(Example):
    is_multiplicative = True
    def __init__(self):
        pass

    def __call__(self, freqs, ants):
        return {ant: i * np.ones_like(freqs) for i, ant in enumerate(ants)}

mult_model = TestMult()
sim.add(mult_model)  # You can also use an instance of the class
np.all(sim.get_data(1,2) == 2)

[47]:

True

[48]:

# You can also give the class an alias it may be referenced by
class TestAlias(Example):
    is_multiplicative = False
    _alias = ("bias",)
    def __init__(self):
        pass

    def __call__(self, lsts, freqs, bias=10):
        return bias * np.ones((lsts.size, freqs.size), dtype=complex)

sim.refresh()
sim.add("bias")  # Or you can just use the name (or an alias) of the model
np.all(sim.data.data_array == 10)

[48]:

True

How does this work? Basically all that happens is that the @component decorator ensures that every model that is subclassed from the new simulation component is tracked in the component’s _models attribute.

[49]:

Example._models

[49]:

{'testadd': __main__.TestAdd,
 'testmult': __main__.TestMult,
 'testalias': __main__.TestAlias,
 'bias': __main__.TestAlias}

It is also worth noting that you can make new models from components that already exist in hera_sim! We could, for example, make a new RFI model like so:

[50]:

class NewRFI(hera_sim.rfi.RFI):
    _alias = ("example_rfi",)
    """Doesn't do anything, just an example."""
    def __init__(self, **kwargs):
        super().__init__(**kwargs)
    def __call__(self, lsts, freqs):
        """This class will be added to the list of known components, though!"""
        return np.zeros((lsts.size, freqs.size), dtype=complex)

In case it isn’t clear: defining a new simulation component and making a model of that component makes hera_sim automatically aware of its existence, and so the Simulator is able to find it with ease (provided there isn’t a name conflict). See the last entry of the following text.

[51]:

print(hera_sim.components.list_all_components())

array:
  lineararray
  hexarray
foreground:
  diffuseforeground | diffuse_foreground
  pointsourceforeground | pntsrc_foreground
noise:
  thermalnoise | thermal_noise
rfi:
  stations | rfi_stations
  impulse | rfi_impulse
  scatter | rfi_scatter
  dtv | rfi_dtv
  newrfi | example_rfi
gain:
  bandpass | gains | bandpass_gain
  reflections | reflection_gains | sigchain_reflections
crosstalk:
  crosscouplingcrosstalk | cross_coupling_xtalk
  crosscouplingspectrum | cross_coupling_spectrum | xtalk_spectrum
  whitenoisecrosstalk | whitenoise_xtalk | white_noise_xtalk
eor:
  noiselikeeor | noiselike_eor
example:
  testadd
  testmult
  testalias | bias

Saving Data in Chunks¶

We finally get to the last of the advanced features of the Simulator: writing data to disk in a way that resembles how the correlator writes HERA data to disk. HERA data files typically only contain a few integrations, and follow a standard naming convention. With the chunk_sim_and_save method, you can write files to disk so that they are chunked into a set number of integrations per file and follow a particular naming scheme. Here’s an example:

[52]:

# We'll use a new temporary directory
tempdir = Path(tempfile.mkdtemp())
sim.chunk_sim_and_save(
    save_dir=tempdir,  # Write the files to tempdir
    Nint_per_file=2,  # Include 2 integrations per file
    prefix="example",  # Prefix the file basename with "example"
    sky_cmp="custom",  # Tack on "custom" after the JD in the filename
    state="true",  # Tack on "true" after "custom" in the filename
    filetype="uvh5",  # Use the uvh5 file format
)

[53]:

print("\n".join(str(f) for f in tempdir.iterdir()))

/tmp/tmpbyom8sqz/example.2458119.50099.custom.true.uvh5
/tmp/tmpbyom8sqz/example.2458119.50173.custom.true.uvh5
/tmp/tmpbyom8sqz/example.2458119.50223.custom.true.uvh5
/tmp/tmpbyom8sqz/example.2458119.50198.custom.true.uvh5
/tmp/tmpbyom8sqz/example.2458119.50124.custom.true.uvh5
/tmp/tmpbyom8sqz/example.2458119.50050.custom.true.uvh5
/tmp/tmpbyom8sqz/example.2458119.50025.custom.true.uvh5
/tmp/tmpbyom8sqz/example.2458119.50000.custom.true.uvh5
/tmp/tmpbyom8sqz/example.2458119.50074.custom.true.uvh5
/tmp/tmpbyom8sqz/example.2458119.50149.custom.true.uvh5

The filename format is save_dir/[{prefix}.]{jd_major}.{jd_minor}[.{sky_cmp}][.{state}].{filetype}. Note that this hasn’t been tested with other filetypes (e.g. miriad), so results may vary if you deviate from the uvh5 format.

It’s also possible to provide reference files for deciding how to chunk the files:

[54]:

sim.chunk_sim_and_save(
    save_dir=tempdir,
    ref_files=list(tempdir.iterdir()),
    prefix="new",
    sky_cmp="example",
    state="files",
    filetype="uvh5",
)

Telescope hera_sim is not in known_telescopes.
Telescope hera_sim is not in known_telescopes.
Telescope hera_sim is not in known_telescopes.
Telescope hera_sim is not in known_telescopes.
Telescope hera_sim is not in known_telescopes.
Telescope hera_sim is not in known_telescopes.
Telescope hera_sim is not in known_telescopes.
Telescope hera_sim is not in known_telescopes.
Telescope hera_sim is not in known_telescopes.
Telescope hera_sim is not in known_telescopes.

[55]:

print("\n".join(str(f) for f in tempdir.iterdir() if str(f).endswith("files.uvh5")))

/tmp/tmpbyom8sqz/new.2458119.50025.example.files.uvh5
/tmp/tmpbyom8sqz/new.2458119.50099.example.files.uvh5
/tmp/tmpbyom8sqz/new.2458119.50149.example.files.uvh5
/tmp/tmpbyom8sqz/new.2458119.50000.example.files.uvh5
/tmp/tmpbyom8sqz/new.2458119.50223.example.files.uvh5
/tmp/tmpbyom8sqz/new.2458119.50074.example.files.uvh5
/tmp/tmpbyom8sqz/new.2458119.50198.example.files.uvh5
/tmp/tmpbyom8sqz/new.2458119.50050.example.files.uvh5
/tmp/tmpbyom8sqz/new.2458119.50124.example.files.uvh5
/tmp/tmpbyom8sqz/new.2458119.50173.example.files.uvh5

The `run_sim` Method¶

The Simulator class also features the run_sim method, which allows you to run an entire simulation with a single method call. The idea behind this is that one might want to decide all of the simulation parameters beforehand, or have a configuration file specifying the simulation parameters, and then run the entire simulation in one go. Below are some examples of how to use the run_sim method.

Defining A Configuration Dictionary¶

We can specify a sequence of steps to be simulated by using a dictionary that maps models to dictionaries that specify their parameter values. This simulation will include diffuse foregrounds, a noiselike EoR signal, and bandpass gains.

[56]:

hera_sim.defaults.deactivate()
sim = Simulator(**sim_params)
config = {
    "diffuse_foreground":
        {
            "Tsky_mdl": hera_sim.defaults("Tsky_mdl"),
            "omega_p": np.ones_like(sim.freqs),
        },
    "noiselike_eor": {"eor_amp": 1e-4},
    "gains": {"dly_rng": (-50, 50)},
}
sim.run_sim(**config)

You have not specified how to seed the random state. This effect might not be exactly recoverable.

[57]:

# Let's take a look at some of the data
fig = waterfall(sim)

FixedFormatter should only be used together with FixedLocator

_images/tutorials_hera_sim_simulator_112_1.png

[58]:

# Now let's verify that the simulation contains what we want.
print(sim.data.history)

hera_sim v1.0.1.dev14+ga590958: Added diffuseforeground using parameters:
Tsky_mdl = <hera_sim.interpolators.Tsky object at 0x7f124cc5b850>
omega_p = [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1.]
hera_sim v1.0.1.dev14+ga590958: Added noiselikeeor using parameters:
eor_amp = 0.0001
hera_sim v1.0.1.dev14+ga590958: Added bandpass using parameters:
dly_rng = (-50, 50)

Using A Configuration File¶

Instead of using a dictionary to specify the parameters, we can instead use a configuration YAML file:

[59]:

cfg_file = tempdir / "config.yaml"
cfg_file.touch()
with open(cfg_file, "w") as cfg:
    cfg.write(
        """
        pntsrc_foreground:
            nsrcs: 10000
            Smin: 0.2
            Smax: 20
            seed: once
        noiselike_eor:
            eor_amp: 0.005
            seed: redundant
        reflections:
            amp: 0.01
            dly: 200
            seed: once
        """
    )
sim.refresh()
sim.run_sim(sim_file=cfg_file)

[60]:

# Again, let's check out some data.
fig = waterfall(sim)

FixedFormatter should only be used together with FixedLocator

_images/tutorials_hera_sim_simulator_117_1.png

[61]:

# Then verify the history.
print(sim.data.history)

hera_sim v1.0.1.dev14+ga590958: Added pointsourceforeground using parameters:
nsrcs = 10000
Smin = 0.2
Smax = 20
hera_sim v1.0.1.dev14+ga590958: Added noiselikeeor using parameters:
eor_amp = 0.005
hera_sim v1.0.1.dev14+ga590958: Added reflections using parameters:
amp = 0.01
dly = 200

[ ]:

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.

[1]:

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:

[2]:

# 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()

[3]:

# 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}

[4]:

# 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.

[5]:

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

[5]:

False

[6]:

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

[6]:

True

[7]:

# 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]

[8]:

# 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

[9]:

# 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()

_images/tutorials_hera_sim_defaults_13_0.png

[10]:

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

[10]:

{'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}}}

[11]:

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

[11]:

{'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}

[12]:

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

[13]:

# 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)

[14]:

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

_images/tutorials_hera_sim_defaults_18_0.png

[15]:

waterfall(sim2, antpairpol)

_images/tutorials_hera_sim_defaults_19_0.png

[ ]:

The following tutorial will give you an overview of how to use hera_sim from the command line to add instrumental systematics / noise to visibilities:

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 (this can be run from anywhere if hera_sim is installed):

$ hera-sim-simulate.py run --help
usage: hera-sim-simulate.py [-h] [-o OUTFILE] [-v] [-sa] [--clobber] config

Run a hera_sim-managed simulation from the command line.

positional arguments:
  config                Path to configuration file.

options:
  -h, --help            show this help message and exit
  -o OUTFILE, --outfile OUTFILE
                        Where to save simulated data. Overrides outfile
                        specified in config.
  -v, --verbose         Print progress updates.
  -sa, --save_all       Save each simulation component.
  --clobber             Overwrite existing files in case of name conflicts.

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
# This document is intended to serve as a template for constructing new
# configuration YAMLs for use with the command-line interface.

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'
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:
   n_freq: 100
   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
   # configuration file or one of the season default keywords
   'h2c'
systematics:
   rfi:
       # 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
           # 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
           Smax: 300
           beta: -1.5
           spectral_index_mean: -1.0
           spectral_index_std: 0.5
           reference_freq: 0.5
       # Note regarding seed_redundantly:
       # This ensures that baselines within a redundant group see the same sky;
       # however, this does not ensure that the sky is actually consistent. So,
       # while the data produced can be absolutely calibrated, it cannot be
       # used to make sensible images (at least, I don't *think* it can be).

simulation:
   # 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,
                rfi,
                sigchain, ]
   # list particular model components to exclude from simulation
   exclude: [sigchain_reflections,
             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.

The following tutorials will give you an overview of how to simulate visibilities from sky models using hera_sim (first from an interpreter, then from the command line):

Visibility Simulator Examples¶

Although hera_sim is primarily aimed at simulating instrumental effects on top of existing visibility data, the package also has the visibilities module that offers a uniform visibility simulation interface to wrappers for several visibility simulators, as well as a few analytic beam models of the HERA antenna.

Visibility Simulation Interface¶

Starting from versions of hera_sim >= 2.0.0, wrappers for VisCPU, healvis, and pyuvsim visibility simulators are provided through the VisCPU, Healvis, and UVSim classes.

The new ModelData object serves as a container for the visiblity data and all information required to perform visibility simulation. A preferred method to initialize this object is to use the classmethod from_config, providing pyuvsim configuration files that specify the observation and telescope parameters, beam model, and sky model. Direct initialization is also supported and currently required to utilize analytic HERA beam models in the hera_sim.beams module (see Manual Initialization). Under the hood, ModelData wraps three objects used for facilitating simulation in pyuvsim, each encapsulating different information required for simulating visibility. * pyuvdata.UVData * A data container for observation parameters, array layout, telescope location, and visibility. * Accesible from the uvdata property. * Package reference: https://pyuvdata.readthedocs.io/en/latest/index.html * pyradiosky.SkyModel * A data container for sky models, including catalog of sources and diffuse sky maps in HEALPix. * Provide methods for conversion from a catalog of point sources to a HEALPix map and vice versa, which is used by the visibility simulator wrappers. * Accesible from the sky_model property. * Package reference: https://pyradiosky.readthedocs.io/en/latest/index.html * pyuvsim.BeamList and the associated beam_ids * Specify the antenna beam models and which beam model each antenna uses. (See Manual Initialization for supported beam types.) * Accesible from the beams and beam_ids properties. * Reference: https://pyuvsim.readthedocs.io/en/latest/classes.html#pyuvsim.BeamList

The VisibilitySimulation class is the new uniform visibility simulation interface. It takes an instance of a VisibilitySimulator class (of which VisCPU, HealVis and UVSim are examples) and a ModelData object as inputs. Visibility simulation can then be executed simply by calling the simulate method on an instance of a VisibilitySimulation class. This returns a Numpy array of visibilities in the same shape as in the UVData.data_array specification, which is also added to the object property VisibilitySimulation.data_model.uvdata.data_array.

The motivation behind the VisibilitySimulation-VisibilitySimulator-ModelData setup is to provide a single uniform interface to many simulators, and the ability to create wrappers for simulators that are not already provided. This interface is heavily based on the pyuvsim configuration specifications and infrastructures, and thus it provides the ability to perform simulation with any simulator given the same static configuration.

Example 1: Demonstrate the Uniform Visibility Simulator Interface¶

[1]:

"""This example simulates visiblities with VisCPU at 10 frequencies and 60
times with the 50 brightest GLEAM sources as a sky model and the HERA Phase I
array as the instrument."""
from hera_sim import DATA_PATH
from hera_sim.visibilities import VisCPU, ModelData, VisibilitySimulation


# Path to the example pyuvsim configuration files.
# These can be found in `hera_sim.DATA_PATH` to follow along.
config_file = (
    DATA_PATH /
    'tutorials_data/visibility_simulator/obsparam_hera_phase1_gleam_top50.yaml'
).as_posix()

# Initialize a ModelData object.
# The pyuvsim configuration files consist of a few yaml and csv files, and a
# sky model file readable by `pyradiosky`. `ModelData.from_config` takes the
# outermost file obsparam_*.yaml as an input.
data_model = ModelData.from_config(config_file)

# Initialize a VisCPU simulator.
simulator = VisCPU()
# simulator = VisCPU(use_pixel_beams=False)
# simulator = UVSim(quiet=True)

# Construct a VisibilitySimulation object.
simulation = VisibilitySimulation(data_model=data_model, simulator=simulator)

# Executing the simulation by calling `simulate` method.
visibility_array = simulation.simulate()

Accessing Simulation Parameters and Visibility¶

All simulation parameters and visibility data are contained in the ModelData object and can be retrieved from the corresponding object attributes and properties.

Example 2: Plot the array layout of the simulation in Example 1¶

[2]:

import matplotlib.pyplot as plt
from pyuvdata import utils as uvutils

uvd = simulation.data_model.uvdata

# Get antennas positions and corresponding antenna numbers
antpos, antnum = uvd.get_ENU_antpos()

# Plot the EN antenna position.
plt.scatter(antpos[:, 0], antpos[:, 1])
for i, antnum in enumerate(uvd.antenna_numbers):
    plt.text(antpos[i, 0], antpos[i, 1], antnum)
plt.show()

_images/tutorials_visibility_simulator_6_0.png

Example 3: Plot the source positions in the SkyModel used in Example 1¶

[3]:

sky_model = simulation.data_model.sky_model
print('Type of sky model:', sky_model.component_type)
print('Number of sources:', sky_model.Ncomponents)

# Extract the source positions and fluxes
ra = sky_model.ra
dec = sky_model.dec
flux = sky_model.stokes[0, 0, :]

# Plot the source positions with color and size showing source fluxes.
plt.scatter(ra, dec, s=flux, c=flux)
plt.colorbar(label='Flux Density [Jy]')
plt.xlabel('RA [degree]')
plt.ylabel('Dec [degree]')
plt.title('50 Brightest GLEAM Sources')
plt.show()

Type of sky model: point
Number of sources: 50

_images/tutorials_visibility_simulator_8_1.png

Example 4: Plot the waterfall of the visibility simulation in Example 1¶

[4]:

import numpy as np

bls = (0, 1, 'xx')
waterfall = uvd.get_data(bls)
extent = [uvd.freq_array[0, 0] / 1e6, uvd.freq_array[0, -1] / 1e6,
          uvd.time_array[-1] - 2458101, uvd.time_array[0] - 2458101]
plt.imshow(np.abs(waterfall), aspect='auto', interpolation='none',
           extent=extent)
plt.xlabel('Frequency [MHz]')
plt.ylabel('JD - 2458101')
plt.colorbar(label=f'Visibility Amplitude [{uvd.vis_units}]')
plt.show()

_images/tutorials_visibility_simulator_10_0.png

For more comprehensive visualization of visibility data (e.g. waterfalls of differences between visibilities, Fourier Transform of a waterfall, and etc), users may be interested in plotting routines in the HERA uvtools package.

Initializing Model Data¶

As mentioned, the ModelData object contains all the data and information required to perform visibility simulation with hera_sim.

The prefered method to initlize a ModelData object is to generate a set pyuvsim configuration files (see https://pyuvsim.readthedocs.io/en/latest/parameter_files.html) and then pass the top-level obsparams_*.yaml file to the classmethod from_config.

data_model = ModelData.from_config(config_file)

If not using configuration files, the underlying objects can be first created and passed to the ModelData constructor for manual initialization.

Users may consider using hera_sim.io.empty_uvdata, pyuvsim.simsetup.initialize_uvdata_from_params, or pyuvsim.simsetup.initialize_uvdata_from_keywords to construct a UVData object.

See the next section on how to construct a sky model with pyradiosky.SkyModel.

The beam models given to the beams parameter can be: i) a pyuvsim.BeamList object, ii) a pyuvdata.UVBeam, or iii) a list of pyuvsim.AnalyticBeam, including its subclasses.

The hera_sim.beams module provides analytic models of the HERA antenna through the PolyBeam, PerturbedPolayBeam, and ZernikeBeam objects (see HERA Memo #081 and #101, and Choudhuri et al 2021). Although these beam models are subclasses of pyuvsim.AnalyticBeam, they are currently not regcogniziable by the pyuvsim configuration file parser, and thus can only be used by initializing ModelData manually.

Example 5: Manaully initialize a `ModelData` object with `PolyBeam`¶

[5]:

from pyuvsim.simsetup import initialize_uvdata_from_params, _complete_uvdata
from pyuvsim import BeamList
from pyradiosky import SkyModel
from hera_sim.beams import PolyBeam

# Initialize just the UVData and beam_ids from the pyuvsim configuraton files.
# The configuration file (same as in Example 1) uses the "airy" type
# `pyuvsim.AnalyticBeam` as the beam models. Here,
# `pyuvsim.simsetup.initialize_uvdata_from_param` returns `UVData`, `BeamList`,
# and `beam_ids`, but we will discard and replace the `BeamList` with a new
# `BeamList` that uses the `PolyBeam` model.
uvdata, _, beam_ids = initialize_uvdata_from_params(config_file)
# Fill `data_array` with zeros and commplete the UVData object.
_complete_uvdata(uvdata, inplace=True)

# Construct a new `BeamList` with polarized `PolyBeam`.
# `beam_coeffs` from Choudhuri et al 2020.
cfg_pol_beam = dict(
    ref_freq=1e8,
    spectral_index=-0.6975,
    beam_coeffs=[
        2.35088101e-01,
        -4.20162599e-01,
        2.99189140e-01,
        -1.54189057e-01,
        3.38651457e-02,
        3.46936067e-02,
        -4.98838130e-02,
        3.23054464e-02,
        -7.56006552e-03,
        -7.24620596e-03,
        7.99563166e-03,
        -2.78125602e-03,
        -8.19945835e-04,
        1.13791191e-03,
        -1.24301372e-04,
        -3.74808752e-04,
        1.93997376e-04,
        -1.72012040e-05,
    ],
    polarized=True,
)
beams = BeamList([PolyBeam(**cfg_pol_beam)])

# Load the SkyModel. Same model used in Example 1.
sky_model_file = (
    DATA_PATH /
    'tutorials_data/visibility_simulator/gleam_top50.skyh5'
).as_posix()
sky_model = SkyModel()
sky_model.read_skyh5(sky_model_file)

# Construct ModelData
data_model_polybeam = ModelData(
    uvdata=uvdata,
    sky_model=sky_model,
    beam_ids=beam_ids,
    beams = beams
)

Constructing a SkyModel¶

pyradiosky.SkyModel is both a container and an interface for representing astrophysical radio sources.

A SkyModel object will usually consist of many components that represent either compact point sources or pixels of a HEALPix map. This choice can be set by specifying component_type="point" or component_type="healpix" when constructing the object. Doing so will determine which parameters are required as followed. * If component_type="point", the name of each source and the source position (ra and dec) are required. * If component_type="healpix", the nside parameter of the HEALPix map, the indices of the HEALPix pixel hpx_inds, and the HEALPix pixel ordering hpx_order must be set. * If component_type is not set, the type is infered from weather nside is given.

The component flux is given per frequency as an array of Stoke IQUV via the stokes argument (Shape: (4, Nfreqs, Ncomponents)). In addition, the spectral_type must be specified, which can be: * “flat” : Flat spectrum. * “full” : Flux is defined at each frequency given in the freq_array attribute. * “subband” : Flux is defined at a set of band centers give in the freq_array attribute. * “spectral_index” : Flux is defined at a reference frequency set in the reference_frequency attribute and follows a power law given a spectral index set in the spectral_index attribute.

The method at_frequencies can be used to evaluate the stokes array of a SkyModel object at specified frequencies, converting the object to a “full” spectral type when inplace=True, or yielding a new instance of the object with a “full” spectral type if inplace=False.

The classmethods pyradiosky.from_gleam_catalog, pyradiosky.from_votable_catalog, and pyradiosky.from_text_catalog can be used to SkyModel object from the GLEAM EGC catalog file, a votable catalog file, or from a tab separated value file. The pyuvsim.simsetup.create_mock_catalog method can also be used to create a mock catalog of point source and/or diffuse sky

Loading a HEALPix map into a SkyModel object is currently only possible via the native skyh5 file format. Thus, a diffuse sky map in a HEALPix FITS format (e.g., Haslam 408 MHz) must be manually initialized from the HEALPix array values and indices.

A method for converting a SkyModel object from one component type to another is provided, which can be used to convert a “healpix” model to “point” and combine with another “point”-type model, for example. This conversion is used by hera_sim in its VisCPU and healvis wrappers.

Below, we show how to create a random point source model and a uniformly distrubuted diffuse model, and combine them into a single SkyModel object. Please refer to the pyradiosky developer API for details of the objects parameters and methods.

Example 6: Construct a random point source SkyModel¶

[6]:

from pyradiosky import SkyModel
from astropy.coordinates import Longitude, Latitude
from astropy import units as u

## Point Soruce Model ##
component_type = 'point'

# Each source is a component in a SkyModel.
# Create a sky model with 50 random sources (i.e. 50 point source components).
nsources = 100

# `ra` and `dec` must be Longitude and Latitude astropy quantity
ra = Longitude(np.random.uniform(0, 360, nsources) * u.degree)
dec = Latitude(np.random.uniform(-90, 90, nsources) * u.degree)

# `names` are required for point soruces
names = [f'rand{i:02d}' for i in range(nsources)]

# SkyModel supports several spectral type. We will use "spectral_index" here,
# where the flux of each source is describe by a power law given a reference
# frequency and a spectral index.
spectral_type = 'spectral_index'
spectral_index = np.random.uniform(-4, 3, nsources)
reference_frequency = 200e6 * np.ones(nsources) * u.Hz

# Source fluxes must be given in Stokes IQUV as an array of astropy quantity
# with an appropriate unit. Shape (4, Nfreqs, Ncomponents).
# Nfreqs is 1 here because we are using "spctral_index".
stokes = np.zeros((4, 1, nsources)) * u.Jy
stokes[0, :, :] = np.random.uniform(10, 100, nsources) * u.Jy
# Let's also give our sources ~10% linar polarization
stokes[1, :, :] = np.random.uniform(1, 5, nsources) * u.Jy
stokes[2, :, :] = np.random.uniform(1, 5, nsources) * u.Jy

# Collect all parameters in a dictionary to be passed to a SkyModel constructor
source_params = {
    'component_type': component_type,
    'name': names,
    'ra': ra,
    'dec': dec,
    'stokes': stokes,
    'spectral_type': spectral_type,
    'reference_frequency': reference_frequency,
    'spectral_index': spectral_index,
    'history': ' Generate a random polarized source catalog with 50 sources.\n'
}

source_model = SkyModel(**source_params)

print('SkyModel History:\n', source_model.history)
print('Type of sky model:', source_model.component_type)
print('Number of components = number of sources:', source_model.Ncomponents)

SkyModel History:
  Generate a random polarized source catalog with 50 sources.
  Read/written with pyradiosky version: 0.1.2.
Type of sky model: point
Number of components = number of sources: 100

Example 7: Construct a uniformly distributed diffuse HEALPix SkyModel¶

[7]:

# Create a diffuse HEALPix SkyModel with a "rainbow" sky and flat spectra.
component_type = 'healpix'
nside = 2 ** 3
npix = 12 * nside ** 2
hpx_order = 'ring'
hpx_inds = np.arange(npix)

spectral_type = 'flat'

stokes = np.zeros((4, 1, npix)) * u.K
stokes[0, :, :] = np.random.uniform(400, 500, npix) * u.K

diffuse_params = {
    'component_type': component_type,
    'nside': nside,
    'hpx_inds': hpx_inds,
    'hpx_order': hpx_order,
    'spectral_type' : spectral_type,
    'stokes': stokes,
    'history': ' Create a diffuse "rainbow" sky with a flat spectra.\n'
}
diffuse_model = SkyModel(**diffuse_params)
print('SkyModel History:\n', diffuse_model.history)
print('Type of sky model:', diffuse_model.component_type)
print('Nside:', diffuse_model.nside)
print('Number of components = number of healpix pixels:',
      diffuse_model.Ncomponents)

SkyModel History:
  Create a diffuse "rainbow" sky with a flat spectra.
  Read/written with pyradiosky version: 0.1.2.
Type of sky model: healpix
Nside: 8
Number of components = number of healpix pixels: 768

Example 8: Convert one SkyModel type to another and concatnate two SkyModel objects¶

One model type can be converted to another type and concatenated to form a composite model. The concatnation can only be done if component_type, spectral_type, and freq_array of both objects matches

[8]:

# Evaluate both sky models at frequencies in Example 1 simulation before
# concatenating.
frequencies = simulation.data_model.uvdata.freq_array[0] * u.Hz
source_model_full = source_model.at_frequencies(frequencies, inplace=False)
diffuse_model_full = diffuse_model.at_frequencies(frequencies, inplace=False)

# Drop spectral_index from the source model before combining.
# See GitHub issue https://github.com/RadioAstronomySoftwareGroup/pyradiosky/issues/160
source_model_full.spectral_index = None

# Convert the diffuse model to "point" type, making sure the flux unit
# is also converted to jansky
diffuse_model_full.healpix_to_point(to_jy=True)
composite_model = source_model_full.concat(diffuse_model_full, inplace=False)

Let’s run some simulations with these sky models and other speciication the same as in Example 1 and compare the output visibilities.

[9]:

# Define a function to run the simulation.
def run_simulation(sky_model):
    """Create and run a simulation given a sky model."""
    # Use the same config file as in Example 1 to initialize UVdata and beams.
    uvdata, beams, beam_ids = initialize_uvdata_from_params(config_file)
    _complete_uvdata(uvdata, inplace=True)

    # Initialize ModelData given the SkyModel
    data_model = ModelData(uvdata=uvdata, sky_model=sky_model,
                           beams=beams, beam_ids=beam_ids)

    # Initialize a VisibilitySimulation instance, using VisCPU
    simulation = VisibilitySimulation(
        data_model=data_model,
        simulator=VisCPU(use_pixel_beams=True)
    )
    _  = simulation.simulate()

    return simulation


source_sim = run_simulation(source_model_full)
diffuse_sim = run_simulation(diffuse_model_full)
composite_sim = run_simulation(composite_model)

[10]:

# Define some labels for all the cases that we want to compare
waterfall_labels = ['A: Source SkyModel', 'B: Diffuse SkyModel',
                    'C: Composite SkyModel',
                    'A + B', 'C - (A + B)']

# Get waterfalls of the complex visibility,
# same baseline (0, 1, 'xx') as in Example 1.
waterfall_list = [sim.data_model.uvdata.get_data(bls)
                  for sim in [source_sim, diffuse_sim, composite_sim]]

# Make A + B, and C - (A + B), add them to the list
waterfall_list.append(waterfall_list[0] + waterfall_list[1])
waterfall_list.append(waterfall_list[2] - waterfall_list[3])

# Extent for the waterfall plot, same for all three simulations
extent = [source_sim.uvdata.freq_array[0, 0] / 1e6,
          source_sim.uvdata.freq_array[0, -1] / 1e6,
          source_sim.uvdata.time_array[-1] - 2458101,
          source_sim.uvdata.time_array[0] - 2458101]

# Plotting
fig, ax = plt.subplots(5, 1, sharex='all', sharey='all', figsize=(6, 8))
for i, waterfall in enumerate(waterfall_list):
    im = ax[i].imshow(np.abs(waterfall), aspect='auto', interpolation='none',
                      extent=extent)
    fig.colorbar(im, ax=ax[i], label='Amplitude [Jy]')
    ax[i].set_title(waterfall_labels[i])
    ax[i].set_ylabel('JD - 2458101')
ax[-1].set_xlabel('Frequency [MHz]')
fig.subplots_adjust(hspace=0.3)

_images/tutorials_visibility_simulator_25_0.png

Simulator Choices¶

hera_sim.visibilities provide wraper classes for three visibility simulators, allowing them to be used in the new uniform interface. These classes are actually subclasses of the VisibilitySimulator class, a new abstract class that serves as a template for wrapping visibility simulators to be used with the uninorm interface, making it easy to create your own wrapper for a different simulator.

Each of the simulator approaches the visibility calculation differently with its own merits.

VisCPU is the wrapper for the vis_cpu package, a Python/numpy-based simulator for interferometer visibilities developed by members from HERA collaboration. It models the sky as an ensemble of point sources, each with their own frequency spectrum. The code is capable of modelling polarized visibilities and primary beams, but currently only a Stokes I sky model.

The HERA Memo #98 describes the mathematical underlying of the vis_cpu simulator although be warned that some information might have already been obsolete due to the recent rapid development of the simulator. In summary, vis_cpu computes a geometric delay for each source and antenna, which is used to calculate the phase factor that is then multiplied by the source flux to obtain per-antenna visibility factor. The direction-dependent polarized antenna beam factor (Jones matrix) is computed by evaluating the complex (efield) antenna pattern at each source position on the azimutal-zenithal coordinates. This is usually done by calling the interp method of the UVBeam object, or by directly evaluating an analytic function in the case of using an analytic beam model. All of these are done per time at a given single frequency.

The hera_sim wrapper is default to used the “pixel beam” mode, which interpolates each anteanna beam once onto an (l, m) grid that is then used to constructs a 2D Bivariate spline. The spline object is evaluated given a source position to estimate the antenna beam pattern. This bypasses the need to compute the new (l, m) coordinates at every time, and, although less accurate, is usually faster for complex efield beams. However, analytic beam models such as the PolyBeam object can precisely and quickly give the beam value at any given az-za coordinate. Thus, pixel beam mode should not be used when using analytic beam models.

The hera_sim wrapper also adds a frequency loop on top of the internal time loop in vis_cpu codes and interfaces between the ModelData to vis_cpu input parameters, including transformation of a HEALPix SkyModel into a point source model. MPI parallelization over frequencies is also implemented on the hera_sim wrapper although shared memory is currently not supported.

The HealVis class provide a wrapper for the healvis simulator. healvis simulates visibility off of HEALPix shells by directly evaluating the radio interferometry measurement equation (RIME) at the healpix pixels, avoiding artifacts from sky projection (see Section 3. of Lanman et al. 2020 for the underlying mathematical formalism). Due to this unique calculation, a point source catalog must be gridded onto a healpix map to use in healvis. The hera_sim wrapper facilites this conversion through the interface between the HealVis wrapper and the ModelData object.

The UVSim class provides a hera_sim wrapper to pyuvsim, a comprehensive simulation package for radio interferometers that emphasizes accuracy and extensibility over speed. Under the hood, it execute pyuvsim.uvsim.run_uvdata_uvsim given the information in the ModelData object.

The VisibilitySimulator abstract base class defines the required class infrastructure for creating a custom simulator.

[11]:

from hera_sim.visibilities import VisibilitySimulator

VisibilitySimulator??

Init signature: VisibilitySimulator()
Source:
class VisibilitySimulator(metaclass=ABCMeta):
    """Base class for all hera_sim compatible visibility simulators."""

    #: Whether this particular simulator has the ability to simulate point
    #: sources directly.
    point_source_ability = True

    #: Whether this particular simulator has the ability to simulate diffuse
    #: maps directly.
    diffuse_ability = False

    __version__ = "unknown"

    @abstractmethod
    def simulate(self, data_model: ModelData) -> np.ndarray:
        """Simulate the visibilities."""
        pass

    def validate(self, data_model: ModelData):
        """Check that the data model complies with the assumptions of the simulator."""
        pass
File:           ~/src/miniconda3/envs/hera_sim_dev/lib/python3.8/site-packages/hera_sim/visibilities/simulators.py
Type:           ABCMeta
Subclasses:     UVSim, VisCPU, HealVis

A custom simulator can be made by subclassing the VisibilitySimulator class and overriding the simulate method, which must take a ModelData object as an input and must return a Numpy array with the same shap as UVData.data_array. The validate method may be override, and point_source_ability and diffuse_abilityatrributes may be set, as neccessary.

Let’s take a look at the codes of the UVSim class that wraps the pyuvsim simulator as an example. It defines a simulate method that calls pyuvsim.uvsim.run_uvdata_uvsim, passing in a UVData object from the ModelData input, and returning the data_array from the UVData output.

[12]:

from hera_sim.visibilities import UVSim

UVSim??

Init signature: UVSim(quiet: bool = False)
Source:
class UVSim(VisibilitySimulator):
    """A wrapper around the pyuvsim simulator.

    Parameters
    ----------
    quiet
        If True, don't print anything.
    """

    def __init__(self, quiet: bool = False):
        self.quiet = quiet

    def simulate(self, data_model: ModelData):
        """Simulate the visibilities."""
        beam_dict = {
            ant: data_model.beam_ids[str(idx)]
            for idx, ant in enumerate(data_model.uvdata.antenna_names)
        }

        # TODO: this can be removed once
        # https://github.com/RadioAstronomySoftwareGroup/pyuvsim/pull/357
        # is mereged.
        if data_model.sky_model.name is not None:
            data_model.sky_model.name = np.array(data_model.sky_model.name)

        warnings.warn(
            "UVSim requires time-ordered data. Ensuring that order in UVData..."
        )
        data_model.uvdata.reorder_blts("time")

        # The UVData object must have correctly ordered pols.
        # TODO: either remove this when pyuvsim fixes bug with ordering
        # (https://github.com/RadioAstronomySoftwareGroup/pyuvsim/issues/370) or
        # at least check whether reordering is necessary once uvdata has that ability.
        if np.any(data_model.uvdata.polarization_array != np.array([-5, -6, -7, -8])):
            warnings.warn(
                "In UVSim, polarization array must be in AIPS order. Reordering..."
            )
            data_model.uvdata.reorder_pols("AIPS")

        out_uv = pyuvsim.uvsim.run_uvdata_uvsim(
            input_uv=data_model.uvdata,
            beam_list=data_model.beams,
            beam_dict=beam_dict,
            catalog=pyuvsim.simsetup.SkyModelData(data_model.sky_model),
            quiet=self.quiet,
        )
        return out_uv.data_array
File:           ~/src/miniconda3/envs/hera_sim_dev/lib/python3.8/site-packages/hera_sim/visibilities/pyuvsim_wrapper.py
Type:           ABCMeta
Subclasses:

[ ]:

Running `hera-sim-vis` from the command line¶

As of v2.4.0 of hera_sim, we have included a command-line interface for performing visibility simulations using any compatible visiblity simulator. The interface for the script is (this can be run from anywhere if hera_sim is installed):

$ hera-sim-vis.py --help
usage: hera-sim-vis.py [-h] [--object_name OBJECT_NAME] [--compress COMPRESS]
                       [--normalize_beams] [--fix_autos]
                       [--max-auto-imag MAX_AUTO_IMAG] [-d] [--run-auto-check]
                       [--profile] [--profile-funcs PROFILE_FUNCS]
                       [--profile-output PROFILE_OUTPUT]
                       [--log-level {INFO,ERROR,WARNING,CRITICAL,DEBUG}]
                       [--log-width LOG_WIDTH] [--log-plain-tracebacks]
                       [--log-absolute-time] [--log-no-mem]
                       [--log-mem-backend {tracemalloc,psutil}]
                       [--log-show-path]
                       obsparam simulator_config

Run vis_cpu via hera_sim given an obsparam.

positional arguments:
  obsparam              pyuvsim-formatted obsparam file.
  simulator_config      YAML configuration file for the simulator.

options:
  -h, --help            show this help message and exit
  --object_name OBJECT_NAME
                        Set object_name in the UVData
  --compress COMPRESS   Compress by redundancy. A file name to store the
                        cache.
  --normalize_beams     Peak normalize the beams.
  --fix_autos           Check and fix non-real xx/yy autos
  --max-auto-imag MAX_AUTO_IMAG
                        Maximum fraction of imaginary/absolute for autos
                        before raising an error
  -d, --dry-run         If set, create the simulator and data model but don't
                        run simulation.
  --run-auto-check      whether to check autos are real

Options for line-profiling:
  --profile             Line-Profile the script
  --profile-funcs PROFILE_FUNCS
                        List of functions to profile
  --profile-output PROFILE_OUTPUT
                        Output file for profiling info.

Options for logging:
  --log-level {INFO,ERROR,WARNING,CRITICAL,DEBUG}
                        logging level to display.
  --log-width LOG_WIDTH
                        width of logging output
  --log-plain-tracebacks
                        use plain instead of rich tracebacks
  --log-absolute-time   show logger time as absolute instead of relative to
                        start
  --log-no-mem          do not show memory usage
  --log-mem-backend {tracemalloc,psutil}
  --log-show-path       show path of code in log msg

Two main configuration files are required. The first is an “obsparam” file, which should follow the formatting defined by pyuvsim. As described in the visibility simulator tutorial, the hera_sim.visibilities module provides a universal interface between this particular configuration setup and a number of particular simulators.

To specify the simulator to be used, the second configuration file must be provided. An example of this configuration file, defined for the VisCPU simulator, can be found in the repo’s config_examples directory. Here are its contents:

$ cat -n ../config_examples/simulator.yaml
simulator: MatVis
precision: 2
ref_time: mean
correct_source_positions: true

Notice that the file contains a simulator: entry. This must be the name of a simulator derived from the base VisibilitySimulator class provided in hera_sim. Usually, this will be one of the built-in simulators listed in the API reference under “Built-In Simulators”.

However, it may also be a custom simulator, defined outside of hera_sim, so long as it inherits from VisibilitySimulator. To use one of these via command line, you need to provide the dot-path to the class, so for example: my_library.module.MySimulator.

The remaining parameters are passed to the constructor for the given simulator. So, for example, for the VisCPU simulator, all available parameters are listed under the Parameters section of its class API. In general, the class may do some transformation of the YAML inputs before constructing its instance, using the from_yaml_dict() method. Check out the documentation for that method for your particular simulator to check if it requires any transformations (eg. if it required data loaded from a file, it might take the filename from the YAML, and read the file in its from_yaml_dict() method before constructing the object with the full data).

API Reference¶

Modules¶

`hera_sim.vis`	Functions for producing white-noise redundant visibilities.
`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`	Visibility-space foreground models.
`hera_sim.interpolators`	This module provides interfaces to different interpolation classes.
`hera_sim.io`	Methods for input/output of data.
`hera_sim.noise`	Models of thermal noise.
`hera_sim.rfi`	Models of radio frequency interference.
`hera_sim.sigchain`	Models of signal-chain systematics.
`hera_sim.simulate`	Module containing a high-level interface for `hera_sim`.
`hera_sim.utils`	Utility module.
`hera_sim.cli_utils`	Useful helper functions and argparsers for running simulations via CLI.
`hera_sim.components`	A module providing discoverability features for hera_sim.

Visibility Simulators¶

Simulation Framework¶

hera_sim.visibilities.simulators

Module defining a high-level visibility simulator wrapper.

Built-In Simulators¶

hera_sim.visibilities.pyuvsim_wrapper.UVSim([quiet])

A wrapper around the pyuvsim simulator.

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:

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

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

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.
Make a development environment. We highly recommend using conda for this. With conda, just run:
```
conda env create -n hera python=3
pip install -e .[dev]
pre-commit install
```

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

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:

Include passing tests (run pytest)
Update documentation when there’s new API, functionality etc.
Add a note to CHANGELOG.rst about the changes.
Add yourself to AUTHORS.rst.

Developing `hera_sim`¶

All docstrings should be written in Numpy docstring format.

Authors¶

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

Changelog¶

dev¶

Added¶

Classes subclassed from SimulationComponent now have a is_randomized class attribute that informs the Simulator of whether it should provide a BitGenerator to the class when simulating the component. - Classes which use a random component should now have a rng attribute,

which should be treated in the same manner as other model parameters. In other words, random states are now effectively treated as model parameters.
New simulator class FFTVis that uses the fftvis package to simulate visibilities. This is a CPU-based visibility simulator that is faster than MatVis for large, compact arrays.

Changed¶

All random number generation now uses the new numpy API. - Rather than seed the global random state, a new BitGenerator is made

with whatever random seed is desired.
- The Simulator API has remained virtually unchanged, but the internal logic that handles random state management has received a significant update.

Deprecated¶

Support for Python 3.9 has been dropped.

Fixed¶

API calls for pyuvdata v2.4.0.

v4.1.0 [2023.06.26]¶

This release heavily focuses on performance of the visibility simulators, and in particular the VisCPU simulator.

Fixed¶

Passing spline_interp_opts now correctly pipes these options through to the visibility simulators.

Added¶

New _blt_order_kws class-attribute for VisibilitySimulator subclasses, that can be used to create the mock metadata in an order corresponding to that required by the simulator (instead of re-ordering after data creation, which can take some time).
New optional compress_data_model() method on VisibilitySimulator subclasses that allows unnecessary metadata in the UVData object to be dropped before simulation (can be restored afterwards with the associated restore_data_model()). This can reduce peak memory usage.
New check_antenna_conjugation parameter for the VisCPU simulator, to allow turning off checks for antenna conjugation, which takes time and is unnecessary for mock datasets.
Dependency on hera-cli-utils which adds options like --log-level and --profile to hera-sim-vis.py.
Option to use a taper in generating a bandpass.
utils.tanh_window function for generating a two-sided tanh window.
interpolators.Reflection class for building a complex reflection coefficient interpolator from a npz archive.
Reflection coefficient and beam integral npz archives for the phase 1 and phase 4 systems (i.e., dipole feed and Vivaldi feed).

Changed¶

run_check_acceptability is now False by default when constructing simulations from obsparams configuration files, to improve performance.
For VisCPU simulator, we no longer copy the whole data array when simulating, but instead just fill the existing one, to save on peak RAM.
Made VisCPU._reorder_vis() much faster (like 99% time reduction).
The --compress option to hera-sim-vis.py is no longer a boolean flag but takes a file argument. This file will be written as a cache of the baseline-time indices required to keep when compressing by redundancy.

v4.0.0 [2023.05.22]¶

Breaking Changes¶

Removed the HealVis wrapper. Use pyuvsim instead.

Changed¶

Updated package to always use future array shapes for pyuvdata objects (this includes updates to PolyBeam and Simulator objects amongst others).

v3.1.1 [2023.02.23]¶

Changed¶

Coupling matrix calculation in MutualCoupling has been updated to correctly calculate the coupling coefficients from the provided E-field beam.

v3.1.0 [2023.01.17]¶

Added¶

MutualCoupling class that simulates the systematic described in Josaitis et al. 2021.
New class attributes for the SimulationComponent class:
- return_type specifies what type of return value to expect;
- attrs_to_pull specifies which Simulator attributes to use.
Some helper functions for MutualCoupling matrix multiplications.
More attributes from the underlying UVData object exposed to the Simulator.

Changed¶

Simulator._update_args logic has been improved.
Simulator attributes lsts, times, and freqs are no longer cached.

v3.0.0¶

Removed¶

Finally removed ability to set use_pixel_beams and bm_pix on the VisCPU simulator. This was removed in v1.0.0 of vis_cpu.
Official support for py37.

Internals¶

Added isort and pyupgrade pre-commit hooks for cleaner code.

v2.3.4 [2022.06.08]¶

Added¶

NotImplementedError raised when trying to simulate noise using an interpolated sky temperature and phase-wrapped LSTs.
More comparison tests of pyuvsim wrapper.

Fixed¶

Inferred integration time in ThermalNoise when phase-wrapped LSTs are used.
Added **kwargs to PolyBeam.interp method to match UVBeam.
healvis wrapper properly sets cross-pol visibilities to zero.

Changed¶

Temporarily forced all UVData objects in the code to use current array shapes.

v2.3.3 [2022.02.21]¶

Added¶

adjustment.interpolate_to_reference now supports interpolating in time when there is a phase wrap in LST.

Changed¶

Some logical statements in adjustment.interpolate_to_reference were changed to use binary operators on logical arrays instead of e.g. np.logical_or.

v2.3.2 [2022.02.18]¶

Added¶

_extract_kwargs attribute added to the SimulationComponent class. This attribute is used by the Simulator to determine which optional parameters should actually be extracted from the data.
antpair optional parameter added to the ThermalNoise class. This is used to determine whether to simulate noise via the radiometer equation (as is appropriate for a cross-correlation) or to just add a bias from the receiver temperature (which is our proxy for what should happen to an auto-correlation).

Fixed¶

The Simulator class now correctly uses the auto-correlations to simulate noise for the cross-correlations.

v2.3.1 [2022.01.19]¶

Fixed¶

Using the normalize_beams option is now possible with the from_config class method.

v2.3.0 [2022.01.19]¶

Added¶

normalize_beams option in ModelData class. Setting this parameter to True enforces peak-normalization on all of the beams used in the simulation. The default behavior is to not peak-normalize the beams.

v2.2.1 [2022.01.14]¶

Added¶

OverAirCrossCoupling now has a parameter amp_norm. This lets the user decide at what distance from the receiverator the gain of the emitted signal is equal to the base amplitude.

Fixed¶

OverAirCrossCoupling now only simulates the systematic for cross-correlations.
ReflectionSpectrum class had its is_multiplicative attribute set to True.

v2.2.0 [2022.01.13]¶

Added¶

New ReflectionSpectrum class to generate multiple reflections over a specified range of delays/amplitudes.

Fixed¶

Corrected some parameter initializations in sigchain module.

v2.1.0 [2022.01.12]¶

Added¶

New OverAirCrossCoupling class to better model crosstalk in H1C data.

Changed¶

Slightly modified Simulator logic for automatically choosing parameter values. This extends the number of cases the class can handle, but will be changed in a future update.

v2.0.0 [2021.11.16]¶

Added¶

New VisibilitySimulator interface. See the `<https://hera-sim.readthedocs.io/en/latest/tutorials/visibility_simulator.html> Visibility Simulator Tutorial`_ for details. This is a breaking change for usage of the visibility simulators, and includes more robust handling of polarization, fixed ordering of data when put back into the UVData objects, more native support for using pyradiosky to define the sky model, and improved support for vis_cpu.
Interface directly to the pyuvsim simulation engine.
Ability to load tutorial data from the installed package.
New and refactored tests for visibility simulations.

Fixed¶

default feed_array for PolyBeam fixed.

Changed¶

Updated tutorial for the visibility simulator interface (see above link).
vis_cpu made an optional extra
removed the conversions module, which is now in the vis_cpu package.
Can now properly use pyuvdata>=2.2.0.

v1.1.1 [2021.08.21]¶

Added¶

Add a Zernike polynomial beam model.

v1.1.0 [2021.08.04]¶

Added¶

Enable polarization support for vis_cpu (handles polarized primary beams, but only Stokes I sky model so far)
Add a polarized version of the analytic PolyBeam model.

v1.0.2 [2021.07.01]¶

Fixed¶

Bug in retrieval of unique LSTs by Simulator when a blt-order other than time-baseline is used has been fixed. LSTs should now be correctly retrieved.
empty_uvdata() now sets the phase_type attribute to “drift”.

v1.0.1 [2021.06.30]¶

Added¶

Fixed¶

Discrepancy in PointSourceForeground documentation and actual implementation has been resolved. Simulated foregrounds now look reasonable.

Changed¶

The time parameters used for generating an example Simulator instance in the tutorial have been updated to match their description.
Simulator tutorial has been changed slightly to account for the foreground fix.

v1.0.0 [2021.06.16]¶

Added¶

adjustment module from HERA Phase 1 Validation work
- adjust_to_reference()
  
  High-level interface for making one set of data comply with another set of data. This may involve rephasing or interpolating in time and/or interpolating in frequency. In the case of a mismatch between the two array layouts, this algorithm will select a subset of antennas to provide the greatest number of unique baselines that remain in the downselected array.
- All other functions in this module exist only to modularize the above function.
cli_utils module providing utility functions for the CLI simulation script.
components module providing an abstract base class for simulation components.
- Any new simulation components should be subclassed from the SimulationComponent ABC. New simulation components subclassed appropriately are automatically discoverable by the Simulator class. A MWE for subclassing new components is as follows:
  @component class Component: pass class Model(Component): ...
  The Component base class tracks any models subclassed from it and makes it discoverable to the Simulator.
New “season” configuration (called "debug"), intended to be used for debugging the Simulator when making changes that might not be easily tested.
chunk_sim_and_save() function from HERA Phase 1 Validation work
- This function allows the user to write a pyuvdata.UVData object to disk in chunks of some set number of integrations per file (either specified directly, or specified implicitly by providing a list of reference files). This is very useful for taking a large simulation and writing it to disk in a way that mimics how the correlator writes files to disk.
Ability to generate noise visibilities based on autocorrelations from the data. This is achieved by providing a value for the autovis parameter in the thermal_noise function (see ThermalNoise).
The vary_gains_in_time() provides an interface for taking a gain spectrum and applying time variation (linear, sinusoidal, or noiselike) to any of the reflection coefficient parameters (amplitude, phase, or delay).
The CrossCouplingSpectrum provides an interface for generating multiple realizations of the cross-coupling systematic spaced logarithmically in amplitude and linearly in delay. This is ported over from the Validation work.

Fixed¶

The reionization signal produced by eor.noiselike_eor is now guaranteed to be real-valued for autocorrelations (although the statistics of the EoR signal for the autocorrelations still need to be investigated for correctness).

Changed¶

API BREAKING CHANGES
- All functions that take frequencies and LSTs as arguments have had their signatures changed to func(lsts, freqs, *args, **kwargs).
- Functions that employ rough_fringe_filter() or rough_delay_filter() as part of the visibility calculation now have parameters delay_filter_kwargs and/or fringe_filter_kwargs, which are dictionaries that are ultimately passed to the filtering functions. foregrounds.diffuse_foreground and eor.noiselike_eor are both affected by this.
- Some parameters have been renamed to enable simpler handling of package-wide defaults. Parameters that have been changed are:
  filter_type -> delay_filter_type in gen_delay_filter()
  
  filter_type -> fringe_filter_type in gen_fringe_filter()
  
  chance -> impulse_chance in rfi_impulse (see Impulse)
  
  strength -> impulse_strength in rfi_impulse (see Impulse)
  
  Similar changes were made in rfi_dtv (DTV) and rfi_scatter (Scatter).
- Any occurrence of the parameter fqs has been replaced with freqs.
- The noise.jy2T function was moved to utils and renamed. See jansky_to_kelvin().
- The parameter fq0 has been renamed to f0 in RfiStation.
- The _listify function has been moved from rfi to utils.
- sigchain.HERA_NRAO_BANDPASS no longer exists in the code, but may be loaded from the file HERA_H1C_BANDPASS.npy in the data directory.
Other Changes
- The Simulator has undergone many changes that make the class much easier to use, while also providing a handful of extra features. The new Simulator provides the following features:
  A universal add() method for applying any of the effects implemented in hera_sim, as well as any custom effects defined by the user.
  
  A get() method that retrieves any previously simulated effect.
  
  The option to apply a simulated effect to only a subset of antennas, baselines, and/or polarizations, accessed through using the vis_filter parameter.
  
  Multiple modes of seeding the random state to achieve a higher degree of realism than previously available.
  
  The calculate_filters() method pre-calculates the fringe-rate and delay filters for the entire array and caches the result. This provides a marginal-to-modest speedup for small arrays, but can provide a significant speedup for very large arrays. Benchmarking results TBD.
  
  An instance of the Simulator may be generated with an empty call to the class if any of the season defaults are active (or if the user has provided some other sufficiently complete set of default settings).
  
  Some of the methods for interacting with the underlying pyuvdata.UVData object have been exposed to the Simulator (e.g. get_data).
  
  An easy reference to the chunk_sim_and_save() function.
- foregrounds, eor, noise, rfi, antpos, and sigchain have been modified to implement the features using callable classes. The old functions still exist for backwards-compatibility, but moving forward any additions to visibility or systematics simulators should be implemented using callable classes and be appropriately subclassed from SimulationComponent.
- empty_uvdata() has had almost all of its parameter values set to default as None. Additionally, the n_freq, n_times, antennas parameters are being deprecated and will be removed in a future release.
- white_noise() is being deprecated. This function has been moved to the utility module and can be found at gen_white_noise().

v0.4.0 [2021.05.01]¶

Added¶

New features added to vis_cpu
- Analytic beam interpolation
  
  Instead of gridding the beam and interpolating the grid using splines, the beam can be interpolated directly by calling its interp method.
  
  The user specifies this by passing use_pixel_beams=False to vis_cpu.
- A simple MPI parallelization scheme
  
  Simulation scripts may be run using mpirun/mpiexec
  
  The user imports mpi4py into their script and passes mpi_comm=MPI.COMM_WORLD to vis_cpu
- New PolyBeam and PerturbedPolyBeam analytic beams (classes)
  
  Derived from pyuvsim.Analytic beam
  
  Based on axisymmetric Chebyshev polynomial fits to the Fagnoni beam.
  
  PerturbedPolyBeam is capable of expressing a range of non-redundancy effects, including per-beam stretch factors, perturbed sidelobes, and ellipticity/rotation.

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

hera_sim¶

Features¶

Documentation¶

Installation¶

Conda users¶

Pip-only install¶

Developer install¶

Versioning¶

Contents¶

Tutorials and FAQs¶

Tour of hera_sim¶

Foregrounds¶

Diffuse Foregrounds¶

Point-Source Foregrounds¶

Diffuse and Point-Source Foregrounds¶

Noise¶

RFI¶

Gains¶

Crosstalk¶

The Simulator Class¶

Setup¶

Introduction¶

Basic Use¶

Overview¶

Adding Effects¶

Retrieving Simulated Effects¶

Saving A Simulation¶

Advanced Use¶

Singling Out Baselines/Antennas/Polarizations¶

Naming Components¶

Pre-Computing Fringe-Rate/Delay Filters¶

Using Custom Models¶

Registering Classes¶

Saving Data in Chunks¶

The run_sim Method¶

Defining A Configuration Dictionary¶

Using A Configuration File¶

Guide for hera_sim Defaults and Simulator¶

Defaults Basics¶

Running hera_sim from the command line¶

BDA¶

Filing¶

Setup¶

Defaults¶

Systematics¶

Sky¶

Simulation¶

Visibility Simulator Examples¶

Visibility Simulation Interface¶

Example 1: Demonstrate the Uniform Visibility Simulator Interface¶

Accessing Simulation Parameters and Visibility¶

Example 2: Plot the array layout of the simulation in Example 1¶

Example 3: Plot the source positions in the SkyModel used in Example 1¶

Example 4: Plot the waterfall of the visibility simulation in Example 1¶

Initializing Model Data¶

Example 5: Manaully initialize a ModelData object with PolyBeam¶

Constructing a SkyModel¶

Example 6: Construct a random point source SkyModel¶

Example 7: Construct a uniformly distributed diffuse HEALPix SkyModel¶

Example 8: Convert one SkyModel type to another and concatnate two SkyModel objects¶

Simulator Choices¶

Running hera-sim-vis from the command line¶

API Reference¶

Modules¶

Visibility Simulators¶

Simulation Framework¶

Built-In Simulators¶

Contributing¶

Bug reports¶

Documentation improvements¶

Feature requests and feedback¶

Development¶

Pull Request Guidelines¶

Developing hera_sim¶

Authors¶

Changelog¶

dev¶

Added¶

Changed¶

Deprecated¶

The `run_sim` Method¶

Running `hera_sim` from the command line¶

Example 5: Manaully initialize a `ModelData` object with `PolyBeam`¶

Running `hera-sim-vis` from the command line¶

Developing `hera_sim`¶