Accessing HERA Phase I Power Spectrum Limit Data

arXiv:2108.02263
Homepage: https://reionization.org/science/public-data-release-1/

Compiled by Nick Kern on behalf of the HERA Collaboration
August, 2021

Description of the data

The data provided are the spherically averaged, dimensionless power spectra $\Delta^2(k)$ for the two frequency subbands and the three fields quoted in HERA et al. 2021. Please see the paper for in-depth details about how the data are reduced and the power spectra are estimated. Included in the files are

• complex-valued power spectra $\Delta^2(k)$ (uvp.data_array)
• thermal noise floor estimate $P_N$ (uvp.stats_array['P_N'])
• power spectra covariance matrix (uvp.cov_array_real)
• power spectra window functions (uvp.window_function_array)

See the associated uvp.get_* methods for accessing the above arrays from a UVPSpec object.

Metadata provided detail things like the amount of time integration, the frequencies used for computing the power spectrum, the cosmology adopted in the normalization, among other things. Below we detail below how to load and plot the data.

%matplotlib inline
# standard imports
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import numpy as np
from scipy import interpolate
import warnings
warnings.filterwarnings('ignore', module='hera_sim')
warnings.filterwarnings('ignore', category=RuntimeWarning)

# HERA-Stack
import hera_pspec as hp
print("hera_pspec version: {}\ngit hash: {}".format(hp.version.version, hp.version.git_hash))

hera_pspec version: 0.3.0
git hash: 0f901a15827af0296eba8e62e40e4ed22f616d9d

def interp_Wcdf(W, k, lower_perc=0.16, upper_perc=0.84):
    """
    Construct CDF from normalized window function and interpolate
    to get k at window func's 16, 50 & 84 percentile.

    Parameters
    ----------
    W : ndarray
        Normalized window function of shape (Nbandpowers, Nk)
    k : ndarray
        vector of k modes of shape (Nk,)

    Returns
    -------
    ndarray
        k of WF's 50th percentile
    ndarray
        dk of WF's 16th (default) percentile from median
    ndarray
        dk of WF's 84th (default) percentile from median
    """
    # get cdf: take sum of only abs(W)
    W = np.abs(W)
    Wcdf = np.array([np.sum(W[:, :i+1].real, axis=1) for i in range(W.shape[1]-1)]).T
    
    # get shifted k such that a symmetric window has 50th perc at max value
    kshift = k[:-1] + np.diff(k) / 2
    
    # interpolate each mode (xvalues are different for each mode!)
    med, low_err, hi_err = [], [], []
    for i, w in enumerate(Wcdf):
        interp = interpolate.interp1d(w, kshift, kind='linear', fill_value='extrapolate')
        m = interp(0.5)  # 50th percentile
        #m = k[np.argmax(W[i])]  # mode
        med.append(m)
        low_err.append(m - interp(lower_perc))
        hi_err.append(interp(upper_perc) - m)

    return np.array(med), np.array(low_err), np.array(hi_err)

# each field is an independent file, containing Band 1 and Band 2
# just load field 1 for now
uvp = hp.UVPSpec()
field = 1
uvp.read_hdf5('pspec_h1c_idr2_field{}.h5'.format(field))

Accessing metadata

There is a lot of metadata attached to the UVPSpec object. Here is an example of how to access some of them. See https://hera-pspec.readthedocs.io/en/latest/ for more info about the UVPSpec object.

# print out general metadata
print("Ntimes:{}, Npols:{}, Nbaseline-pairs:{}, Nspws:{}".format(uvp.Ntimes, uvp.Npols, uvp.Nblpairs, uvp.Nspws))
print("units:{}".format(uvp.units))

Ntimes:1, Npols:1, Nbaseline-pairs:1, Nspws:2
units:(mK)^2 h^-3 Mpc^3 k^3 / (2pi^2)

# getting the cosmology
print(uvp.cosmo)

Cosmo_Conversions object at <0x7ffac873fdf0>
Om_L : 0.6844; Om_b : 0.0491; Om_c : 0.2644; Om_M : 0.3135; Om_k : 0.0021; H0 : 67.2700

# spectral windows: (start_freq Hz, end_freq Hz, N_freq_channels, N_kbins)
spw_ranges = uvp.get_spw_ranges()
print(spw_ranges)

[(117089843.75, 132714843.75, 160, 40), (150292968.75, 167871093.75, 180, 40)]

# spw metadata
for spw in uvp.spw_array:
    print("\nBand {}\n".format(spw+1) + '-'*15)
    avg_f = np.mean(spw_ranges[spw][:2])
    z = uvp.cosmo.f2z(avg_f)
    effective_bandwidth = float(np.diff(spw_ranges[spw][:2])) / 2
    print("band-center z : {:.2f}".format(z))
    print("band-center freq: {:.2f} MHz".format(avg_f/1e6))
    print("effective bandwidth : {:.1f} MHz".format(effective_bandwidth/1e6))

Band 1
---------------
band-center z : 10.37
band-center freq: 124.90 MHz
effective bandwidth : 7.8 MHz

Band 2
---------------
band-center z : 7.93
band-center freq: 159.08 MHz
effective bandwidth : 8.8 MHz

Accessing the data

The data are accessed via a key, formed as (spw, baseline-pair, polarization). For the fully averaged data, there are only two available keys, which corresponds to the two bands associated with this field.

The data, errorbars, and window functions all have uvp.get_ methods that take one of these keys. We will demonstrate below.

# print the two available keys
band1_key, band2_key = uvp.get_all_keys()
keys = [band1_key, band2_key]
print("band1_key : {}\nband2_key : {}".format(band1_key, band2_key))

band1_key : (0, ((1, 14), (12, 26)), ('pI', 'pI'))
band2_key : (1, ((1, 14), (12, 26)), ('pI', 'pI'))

# get data
band1_dsq = uvp.get_data(band1_key)
band1_cov = uvp.get_cov(band1_key)
band1_wfn = uvp.get_window_function(band1_key)
print("data shape : {}, (N_times, N_kbins)".format(band1_dsq.shape))
print("cov shape: {}, (N_times, N_kbins, N_kbins)".format(band1_cov.shape))
print("win func shape : {}, (N_times, N_kbins, N_kbins)".format(band1_wfn.shape))

data shape : (1, 40), (N_times, N_kbins)
cov shape: (1, 40, 40), (N_times, N_kbins, N_kbins)
win func shape : (1, 40, 40), (N_times, N_kbins, N_kbins)

# get data
band2_dsq = uvp.get_data(band2_key)
band2_cov = uvp.get_cov(band2_key)
band2_wfn = uvp.get_window_function(band2_key)
print("data shape : {}, (N_times, N_kbins)".format(band2_dsq.shape))
print("cov shape: {}, (N_times, N_kbins, N_kbins)".format(band2_cov.shape))
print("win func shape : {}, (N_times, N_kbins, N_kbins)".format(band2_wfn.shape))

data shape : (1, 40), (N_times, N_kbins)
cov shape: (1, 40, 40), (N_times, N_kbins, N_kbins)
win func shape : (1, 40, 40), (N_times, N_kbins, N_kbins)

# plot real component of power spectrum
grid = gridspec.GridSpec(3, 3)
fig = plt.figure(figsize=(16, 8))
fig.subplots_adjust(hspace=0)

spw = 0
kbins = uvp.get_kparas(spw)  # after spherical binning, k_perp=0 so k_mag = k_para

ks = slice(2, None)
xlim = (0, 2.0)
band1_err = np.sqrt(band1_cov[0].diagonal())

y = band1_dsq[0, ks].real
ypos = y > 0
yerr = band1_err[ks].real
x, xerr_low, xerr_hi = interp_Wcdf(band1_wfn[0], kbins, lower_perc=0.16, upper_perc=0.84)
xerr = np.array([xerr_low, xerr_hi]).T[ks]

ax = fig.add_subplot(grid[:2, :])
ax.tick_params(labelsize=15, direction='in', size=5)
ax.tick_params(direction='in', size=3, which='minor')
ax.errorbar(kbins[ks][ypos], y[ypos], c='deeppink', ls='', marker='o', ms=10, yerr=2*yerr[ypos], xerr=xerr[ypos].T)
ax.errorbar(kbins[ks][~ypos], np.abs(y[~ypos])*0, c='deeppink', ls='', marker='o', ms=10, mfc='w', yerr=2*yerr[~ypos])

ax.set_xlim(xlim)
ax.set_ylim(2e1, 1e7)
ax.set_yscale('log')
ax.set_xlabel(r'$k\ [h\ {\rm Mpc}^{-1}]$', fontsize=20)
ax.set_ylabel(r'$\Delta^2(k)\ [{\rm mK}^2]$', fontsize=20)
ax.set_xticklabels([])
ax.set_title("Band {}, Field {}".format(spw+1, field))

ax = fig.add_subplot(grid[2:, :])
ax.tick_params(labelsize=15, direction='in', size=5)
ax.plot(kbins[ks], band1_wfn[0, ks, ks], c='k', lw=1.5)
ax.set_xlim(xlim)
ax.set_ylim(0, 0.8)
ax.set_xlabel(r'$k\ [h\ {\rm Mpc}^{-1}]$', fontsize=20)
ax.set_ylabel(r"$W(k)$", fontsize=20);

Figure 1 | Band 1, Field 1 power spectra. Vertical errorbars are $2\sigma$, horizontal errorbars are the 16 and 84 percentile of the window functions. Negative measured bandpowers are not shown.

# plot real component of power spectrum
grid = gridspec.GridSpec(3, 3)
fig = plt.figure(figsize=(16, 8))
fig.subplots_adjust(hspace=0)

spw = 1
kbins = uvp.get_kparas(spw)

ks = slice(2, None)
xlim = (0, 2.0)
band2_err = np.sqrt(band2_cov[0].diagonal())

y = band2_dsq[0, ks].real
ypos = y > 0
yerr = band2_err[ks].real
x, xerr_low, xerr_hi = interp_Wcdf(band2_wfn[0], kbins, lower_perc=0.16, upper_perc=0.84)
xerr = np.array([xerr_low, xerr_hi]).T[ks]

ax = fig.add_subplot(grid[:2, :])
ax.tick_params(labelsize=15, direction='in', size=5)
ax.tick_params(direction='in', size=3, which='minor')
ax.errorbar(kbins[ks][ypos], y[ypos], c='deeppink', ls='', marker='o', ms=10, yerr=2*yerr[ypos], xerr=xerr[ypos].T)
ax.errorbar(kbins[ks][~ypos], np.abs(y[~ypos])*0, c='deeppink', ls='', marker='o', ms=10, mfc='w', yerr=2*yerr[~ypos])

ax.set_xlim(xlim)
ax.set_ylim(2e1, 1e7)
ax.set_yscale('log')
ax.set_xlabel(r'$k\ [h\ {\rm Mpc}^{-1}]$', fontsize=20)
ax.set_ylabel(r'$\Delta^2(k)\ [{\rm mK}^2]$', fontsize=20)
ax.set_xticklabels([])
ax.set_title("Band {}, Field {}".format(spw+1, field))

ax = fig.add_subplot(grid[2:, :])
ax.tick_params(labelsize=15, direction='in', size=5)
ax.plot(kbins[ks], band2_wfn[0, ks, ks], c='k', lw=1.5)
ax.set_xlim(xlim)
ax.set_ylim(0, 0.8)
ax.set_xlabel(r'$k\ [h\ {\rm Mpc}^{-1}]$', fontsize=20)
ax.set_ylabel(r"$W(k)$", fontsize=20);

Figure 2 | Band 2, Field 1 power spectra. Vertical errorbars are $2\sigma$, horizontal errorbars are the 16 and 84 percentile of the window functions. Negative measured bandpowers are not shown.