HERA delay power spectrum estimation

The hera_pspec library provides all of the tools and data structures needed to perform a delay spectrum analysis on interferometric data. The input data can be in any format supported by pyuvdata, and the output data are stored in HDF5 containers.

You can find the code in the hera_pspec GitHub repository. A set of example Jupyter notebooks are also available on GitHub.

Power spectrum calculations

The PSpecData class takes a set of UVData objects containing visibilities and calculates delay power spectra from them. These are output as UVPSpec objects.

Example delay power spectrum calculation

The following example shows how to load UVData objects into a PSpecData object, specify a set of baselines and datasets that should be cross-multiplied, specify a set of spectral windows, weights, and tapers, and output a set of delay power spectra into a UVPSpec object.

# Load into UVData objects
uvd1 = UVData(); uvd2 = UVData()
uvd1.read_miriad(datafile1)
uvd2.read_miriad(datafile2)

# Create a new PSpecData object
ds = hp.PSpecData(dsets=[uvd1, uvd2], beam=beam)

# bls1 and bls2 are lists of tuples specifying baselines (antenna pairs)
# Here, we specify three baseline-pairs, i.e. bls1[0] x bls2[0],
# bls1[1] x bls2[1], and bls1[2] x bls2[2].
bls1 = [(24,25), (37,38), (38,39)]
bls2 = [(37,38), (38,39), (24,25)]

# Calculate cross-spectra of visibilities for baselines bls1[i] x bls2[i],
# where bls1 are the baselines to take from dataset 0 and bls2 are taken from
# dataset 1 (and i goes from 0..2). This is done for two spectral windows
# (freq. channel indexes between 300-400 and 600-721), with unity weights
# and a Blackman-Harris taper in each spectral window
uvp = ds.pspec(bls1, bls2, dsets=(0, 1), spw_ranges=[(300, 400), (600, 721)],
               input_data_weight='identity', norm='I', taper='blackman-harris',
               verbose=True)

uvp is now a UVPSpec object containing 2 x 3 x Ntimes delay spectra, where 3 is the number of baseline-pairs (i.e. len(bls1) == len(bls2) == 3), 2 is the number of spectral windows, and Ntimes is the number of LST bins in the input UVData objects. Each delay spectrum has length Nfreqs, i.e. the number of frequency channels in each spectral window.

To get power spectra from the UVPSpec object that was returned by pspec:

# Key specifying desired spectral window, baseline-pair, and polarization pair
spw = 1
polpair = ('xx', 'xx')
blpair =((24, 25), (37, 38))
key = (spw, blpair, polpair)

# Get delay bins and power spectrum
dlys = uvp.get_dlys(spw)
ps = uvp.get_data(key)

PSpecData: Calculate optimal quadratic estimates of delay power spectra

The PSpecData class implements an optimal quadratic estimator for delay power spectra. It takes as its inputs a set of UVData objects containing visibilities, plus objects containing supporting information such as weights/flags, frequency-frequency covariance matrices, and PSpecBeam: Primary beam models.

Once data have been loaded into a PSpecData object, the pspec() method can be used to calculate delay spectra for any combination of datasets, baselines, spectral windows etc. that you specify. Some parts of the calculation (e.g. empirical covariance matrix estimation) are cached within the PSpecData object to help speed up subsequent calls to pspec().

Note

The input datasets should have compatible shapes, i.e. contain the same number of frequency channels and LST bins. The validate_datasets() method (automatically called by pspec()) checks for compatibility, and will raise warnings or exceptions if problems are found. You can use the pyuvdata.UVData.select() method to select out compatible chunks of UVData files if needed.

Specifying which spectra to calculate

Each call to pspec() must specify a set of baseline-pairs, a set of datasets, and a set of spectral windows that the power spectrum should be estimated for.

• Datasets correspond to the UVData objects that were stored inside the PSpecData object, and are identified either by an index (numbered according to the order that they were added to the PSpecData object), or label strings (if you specified labels when you added the datasets). A pair of datasets is then specified using the dsets argument, e.g. dsets=(0, 1) corresponds to the first and second datasets added to the PSpecData object. You can specify the same dataset if you want, e.g. dsets=(1, 1), although you should beware of noise bias in this case.
• Baseline pairs are specified as two lists: bls1 is the list of baselines from the first dataset in the pair specified by the dsets argument, and bls2 is the list from the second. The baseline pairs are formed by matching each element from the first list with the corresponding element from the second, e.g. blpair[i] = bls1[i] x bls2[i]. A couple of helper functions are provided to construct appropriately paired lists to calculate all of the cross-spectra within a redundant baseline group, for example; see construct_blpairs() and validate_bls().
• Spectral windows are specified as a list of tuples using the spw_ranges argument, with each tuple specifying the beginning and end frequency channel of the spectral window. For example, spw_ranges=[(220, 320)] defines a spectral window from channel 220 to 320, as indexed by the UVData objects. The larger the spectral window, the longer it will take to calculate the power spectra. Note that
• Polarizations are looped over by default. At the moment, pspec() will only compute power spectra for matching polarizations from datasets 1 and 2. If the UVData objects stored inside the PSpecData object have incompatible polarizations, validate_datasets() will raise an exception.

Note

If the input datasets are phased slightly differently (e.g. due to offsets in LST bins), you can rephase (fringe-stop) them to help reduce decoherence by using the rephase_to_dset() method. Note that the validate_datasets() method automatically checks for discrepancies in how the UVData objects are phased, and will raise warnings or errors if any problems are found.

The PSpecData class

The most frequently-used methods from PSpecData are listed below. See PSpecData Class for a full listing of all methods provided by PSpecData.

class hera_pspec.PSpecData(dsets=[], wgts=None, dsets_std=None, labels=None, beam=None, cals=None, cal_flag=True)

__init__(dsets=[], wgts=None, dsets_std=None, labels=None, beam=None, cals=None, cal_flag=True)

Object to store multiple sets of UVData visibilities and perform operations such as power spectrum estimation on them.

Parameters:

• dsets (list or dict of UVData objects, optional) – Set of UVData objects containing the data that will be used to compute the power spectrum. If specified as a dict, the key names will be used to tag each dataset. Default: Empty list.
• dsets_std (list or dict of UVData objects, optional) – Set of UVData objects containing the standard deviations of each data point in UVData objects in dsets. If specified as a dict, the key names will be used to tag each dataset. Default: [].
• wgts (list or dict of UVData objects, optional) – Set of UVData objects containing weights for the input data. Default: None (will use the flags of each input UVData object).
• labels (list of str, optional) – An ordered list of names/labels for each dataset, if dsets was specified as a list. If None, names will not be assigned to the datasets. If dsets was specified as a dict, the keys of that dict will be used instead of this. Default: None.
• beam (PspecBeam object, optional) – PspecBeam object containing information about the primary beam Default: None.
• cals (list of UVCal objects, optional) – Calibration objects to apply to data. One per dset or one for all dsets.
• cal_flag (bool, optional) – If True, propagate flags from calibration into data

add(dsets, wgts, labels=None, dsets_std=None, cals=None, cal_flag=True)

Add a dataset to the collection in this PSpecData object.

Parameters:

• dsets (UVData or list or dict) – UVData object or list of UVData objects containing data to add to the collection.
• wgts (UVData or list or dict) – UVData object or list of UVData objects containing weights to add to the collection. Must be the same length as dsets. If a weight is set to None, the flags of the corresponding dset are used.
• labels (list of str) – An ordered list of names/labels for each dataset, if dsets was specified as a list. If dsets was specified as a dict, the keys of that dict will be used instead.
• dsets_std (UVData or list or dict) – Optional UVData object or list of UVData objects containing the standard deviations (real and imaginary) of data to add to the collection. If dsets is a dict, will assume dsets_std is a dict and if dsets is a list, will assume dsets_std is a list.
• cals (UVCal or list, optional) – UVCal objects to apply to data.
• cal_flag (bool, optional) – If True, propagate flags from calibration into data

delays()

Return an array of delays, tau, corresponding to the bins of the delay power spectrum output by pspec() using self.spw_range to specify the spectral window.

Returns:delays – Delays, tau. Units: ns. Return type:array_like

pspec(bls1, bls2, dsets, pols, n_dlys=None, input_data_weight='identity', norm='I', taper='none', sampling=False, little_h=True, spw_ranges=None, symmetric_taper=True, baseline_tol=1.0, store_cov=False, store_cov_diag=False, return_q=False, store_window=True, verbose=True, filter_extensions=None, exact_norm=False, history='', r_params=None, cov_model='empirical', known_cov=None)

Estimate the delay power spectrum from a pair of datasets contained in this object, using the optimal quadratic estimator of arXiv:1502.06016.

In this formulation, the power spectrum is proportional to the visibility data via

P = M data_{LH} E data_{RH}

where E contains the data weighting and FT matrices, M is a normalization matrix, and the two separate datasets are denoted as “left-hand” and “right-hand”.

Each power spectrum is generated by taking a baseline (specified by an antenna-pair and polarization key) from bls1 out of dsets[0] and assigning it as data_LH, and a bl from bls2 out of the dsets[1] and assigning it as data_RH.

If the bl chosen from bls1 is (ant1, ant2) and the bl chosen from bls2 is (ant3, ant4), the “baseline-pair” describing their cross multiplication is ((ant1, ant2), (ant3, ant4)).

Parameters:

• bls1, bls2 (list) – List of baseline groups, each group being a list of ant-pair tuples.

• dsets (length-2 tuple or list) – Contains indices of self.dsets to use in forming power spectra, where the first index is for the Left-Hand dataset and second index is used for the Right-Hand dataset (see above).

• pols (tuple or list of tuple) – Contains polarization pairs to use in forming power spectra e.g. (‘XX’,’XX’) or [(‘XX’,’XX’),(‘pI’,’pI’)] or a list of polarization pairs. Individual strings are also supported, and will be expanded into a matching pair of polarizations, e.g. ‘xx’ becomes (‘xx’, ‘xx’).

  If a primary_beam is specified, only equal-polarization pairs can be cross-correlated, as the beam scalar normalization is only implemented in this case. To obtain unnormalized spectra for pairs of different polarizations, set the primary_beam to None.

• n_dlys (list of integer, optional) – The number of delay bins to use. The order in the list corresponds to the order in spw_ranges. Default: None, which then sets n_dlys = number of frequencies.

• input_data_weight (str, optional) – String specifying which weighting matrix to apply to the input data. See the options in the set_R() method for details. Default: ‘identity’.

• norm (str, optional) – String specifying how to choose the normalization matrix, M. See the ‘mode’ argument of get_MW() for options. Default: ‘I’.

• taper (str, optional) – Tapering (window) function to apply to the data. Takes the same arguments as uvtools.dspec.gen_window(). Default: ‘none’.

• sampling (boolean, optional) – Whether output pspec values are samples at various delay bins or are integrated bandpowers over delay bins. Default: False

• little_h (boolean, optional) – Whether to have cosmological length units be h^-1 Mpc or Mpc Default: h^-1 Mpc

• spw_ranges (list of tuples, optional) – A list of spectral window channel ranges to select within the total bandwidth of the datasets, each of which forms an independent power spectrum estimate. Example: [(220, 320), (650, 775)]. Each tuple should contain a start (inclusive) and stop (exclusive) channel used to index the freq_array of each dataset. The default (None) is to use the entire band provided in each dataset.

• symmetric_taper (bool, optional) – Specify if taper should be applied symmetrically to K-matrix (if true) or on the left (if False). Default: True

• baseline_tol (float, optional) – Distance tolerance for notion of baseline “redundancy” in meters. Default: 1.0.

• store_cov (bool, optional) – If True, calculate an analytic covariance between bandpowers given an input visibility noise model, and store the output in the UVPSpec object.

• store_cov_diag (bool, optional) – If True, store the square root of the diagonal of the output covariance matrix calculated by using get_analytic_covariance(). The error bars will be stored in the form of: sqrt(diag(cov_array_real)) + 1.j*sqrt(diag(cov_array_imag)). It’s a way to save the disk space since the whole cov_array data with a size of Ndlys x Ndlys x Ntimes x Nblpairs x Nspws is too large.

• return_q (bool, optional) – If True, return the results (delay spectra and covariance matrices) for the unnormalized bandpowers in the UVPSpec object.

• store_window (bool, optional) – If True, store the window function of the bandpowers. Default: True

• cov_model (string, optional) – Type of covariance model to calculate, if not cached.

  Options=[‘empirical’, ‘dsets’, ‘autos’, ‘foreground_dependent’, (other model names in known_cov)]

  In ‘dsets’ mode, error bars are estimated from user-provided per baseline and per channel standard deivations.

  In ‘empirical’ mode, error bars are estimated from the data by averaging the channel-channel covariance of each baseline over time and then applying the appropriate linear transformations to these frequency-domain covariances.

  In ‘autos’ mode, the covariances of the input data over a baseline is estimated from the autocorrelations of the two antennas forming the baseline across channel bandwidth and integration time.

  In ‘foreground_dependent’ mode, it involves using auto-correlation amplitudes to model the input noise covariance and visibility outer products to model the input systematics covariance.

  For more details see ds.get_analytic_covariance().

• known_cov (dicts of input covariance matrices) – known_cov has a type {Ckey:covariance}, which is the same as ds._C. The matrices stored in known_cov are constructed outside the PSpecData object, unlike those in ds._C which are constructed internally.

  The Ckey should conform to: (dset_pair_index, blpair_int, model, time_index, conj_1, conj_2), e.g. ((0, 1), ((25,37,”xx”), (25, 37, “xx”)), ‘empirical’, False, True), while covariance are ndarrays with shape (Nfreqs, Nfreqs).

  Also see PSpecData.set_C() for more details.

• verbose (bool, optional) – If True, print progress, warnings and debugging info to stdout.

• filter_extensions (list of 2-tuple or 2-list, optional) – Set number of channels to extend filtering width.

• exact_norm (bool, optional) – If True, estimates power spectrum using Q instead of Q_alt (HERA memo #44). The default options is False. Beware that computing a power spectrum when exact_norm is set to False runs two times faster than setting it to True.

• history (str, optional) – History string to attach to UVPSpec object

• r_params (dictionary with parameters for weighting matrix.) – Proper fields and formats depend on the mode of data_weighting.

  For data_weighting set to ‘dayenu’, r_params should be a dict with the following fields:

  filter_centers, a list of floats (or float) specifying the (delay) channel numbers at which to center filtering windows. Can specify fractional channel number.

  filter_half_widths, a list of floats (or float) specifying the width of each filter window in (delay) channel numbers. Can specify fractional channel number.

  filter_factors, list of floats (or float) specifying how much power within each filter window is to be suppressed.

  Absence of an r_params dictionary will result in an error.

Returns:

uvp – Instance of UVPSpec that holds the normalized output power spectrum data.

Return type:

UVPSpec object

Examples

Example 1: No grouping; i.e. each baseline is its own group, no brackets needed for each bl. If:

A = (1, 2); B = (2, 3); C = (3, 4); D = (4, 5); E = (5, 6); F = (6, 7)

and:

bls1 = [ A, B, C ]
bls2 = [ D, E, F ]

then:

blpairs = [ (A, D), (B, E), (C, F) ]

Example 2: Grouping; blpairs come in lists of blgroups, which are considered “grouped” in OQE. If:

bls1 = [ [A, B], [C, D] ]
bls2 = [ [C, D], [E, F] ]

then:

blpairs = [ [(A, C), (B, D)], [(C, E), (D, F)] ]

Example 3: Mixed grouping; i.e. some blpairs are grouped, others are not. If:

bls1 = [ [A, B], C ]
bls2 = [ [D, E], F ]

then:

blpairs = [ [(A, D), (B, E)], (C, F)]

rephase_to_dset(dset_index=0, inplace=True)

Rephase visibility data in self.dsets to the LST grid of dset[dset_index] using hera_cal.utils.lst_rephase.

Each integration in all other dsets is phased to the center of the corresponding LST bin (by index) in dset[dset_index].

Will only phase if the dataset’s phase type is ‘drift’. This is because the rephasing algorithm assumes the data is drift-phased when applying the phasor term.

Note that PSpecData.Jy_to_mK() must be run after rephase_to_dset(), if one intends to use the former capability at any point.

Parameters:

• dset_index (int or str) – Index or label of dataset in self.dset to phase other datasets to.
• inplace (bool, optional) – If True, edits data in dsets in-memory. Else, makes a copy of dsets, edits data in the copy and returns to user.

Returns:

• if inplace – return new_dsets
• else – return None

scalar(polpair, little_h=True, num_steps=2000, beam=None, taper_override='no_override', exact_norm=False)

Computes the scalar function to convert a power spectrum estimate in “telescope units” to cosmological units, using self.spw_range to set spectral window.

See arxiv:1304.4991 and HERA memo #27 for details.

This function uses the state of self.taper in constructing scalar. See PSpecData.pspec for details.

Parameters:

• polpair (tuple, int, or str) – Which pair of polarizations to compute the beam scalar for, e.g. (‘pI’, ‘pI’) or (‘XX’, ‘YY’). If string, will assume that the specified polarization is to be cross-correlated with itself, e.g. ‘XX’ implies (‘XX’, ‘XX’).
• little_h (boolean, optional) – Whether to have cosmological length units be h^-1 Mpc or Mpc Default: h^-1 Mpc
• num_steps (int, optional) – Number of steps to use when interpolating primary beams for numerical integral Default: 10000
• beam (PSpecBeam object) – Option to use a manually-fed PSpecBeam object instead of using self.primary_beam.
• taper_override (str, optional) – Option to override the taper chosen in self.taper (does not overwrite self.taper; just applies to this function). Default: no_override
• exact_norm (boolean, optional) – If True, scalar would just be the X2Y term, as the beam and spectral terms are taken into account while constructing power spectrum.

Returns:

scalar – [int dnu (Omega_PP / Omega_P^2) ( B_PP / B_P^2 ) / (X^2 Y)]^-1 in h^-3 Mpc^3 or Mpc^3.

Return type:

float

units(little_h=True)

Return the units of the power spectrum. These are inferred from the units reported by the input visibilities (UVData objects).

Parameters:little_h (boolean, optional) – Whether to have cosmological length units be h^-1 Mpc or Mpc Default: h^-1 Mpc Returns:pspec_units – Units of the power spectrum that is returned by pspec(). Return type:str

UVPSpec: Container for power spectra

The pspec() method outputs power spectra as a single UVPSpec object, which also contains metadata and various methods for accessing the data, input/output etc.

To access the power spectra, use the get_data() method, which takes a key of the form: (spw, blpair, polpair). For example:

# Key specifying desired spectral window, baseline-pair, and polarization
spw = 1
polpair = ('xx', 'xx')
blpair =((24, 25), (37, 38))
key = (spw, blpair, polpair)

# Get delay bins and power spectrum
dlys = uvp.get_dlys(spw)
ps = uvp.get_data(key)

A number of methods are provided for returning useful metadata:

• get_integrations(): Get the average integration time (in seconds) for a given delay spectrum.
• get_nsamples(): If the power spectra have been incoherently averaged (i.e. averaged after squaring), this is the effective number of samples in that average.
• get_dlys(): Get the delay for each bin of the delay power spectra (in seconds).
• get_blpair_seps(): Get the average baseline separation for a baseline pair, in the ENU frame, in meters.

Dimensions and indexing of the UVPSpec data array

The power spectra are stored internally in UVPSpec.data_array, which is a list of three-dimensional numpy arrays, one for each spectral window. Spectral window indices can be retrieved using the spw_to_indices() method. Each 3D array has shape (Nblpairts, Ndlys, Npols).

• Npols is the number of polarizations. Available polarizations can be retrieved from the UVPSpec.pol_array attribute. This dimension can be indexed using the pol_to_indices() method.

• Ndlys is the number of delays, which is equal to the number of frequency channels within the spectral window. The available frequencies/delays can be retrievd from the UVPSpec.freq_array and UVPSpec.dly_array attributes. Alternatively, use the get_dlys() method to get the delay values.

• Nblpairts is the number of unique combinations of baseline-pairs and times (or equivalently LSTs), i.e. the total number of delay spectra that were calculated for a given polarization and spectral window. Baseline-pairs and times have been collapsed into a single dimension because each baseline-pair can have a different number of time samples.

  You can access slices of the baseline-pair/time dimension using the blpair_to_indices() and time_to_indices() methods. The baseline-pairs and times contained in the object can be retrieved from the UVPSpec.blpair_array and UVPSpec.time_avg_array (or UVPSpec.lst_avg_array) attributes.

Note

The UVPSpec.time_avg_array attribute is just the average of the times of the input datasets. To access the original times from each dataset, see the UVPSpec.time_1_array and UVPSpec.time_2_array attributes (or equivalently UVPSpec.lst_1_array and UVPSpec.lst_2_array).

Averaging and folding spectra

By default, separate delay spectra are produced for every LST bin and polarization available in the input datasets, and for every baseline-pair passed to pspec(). To (incoherently) average these down into a single delay spectrum, use the average_spectra() method. For example, to average over all times and baseline-pairs (i.e. assuming that the UVPSpec contains spectra for a single redundant baseline group):

# Build a list of all baseline-pairs in the UVPSpec object
blpair_group = [sorted(np.unique(uvp.blpair_array))]

# Average over all baseline pairs and times, and return in a new ``UVPSpec`` object
uvp_avg = uvp.average_spectra(blpair_groups=blpair_group, time_avg=True, inplace=False)

For UVPSpec objects containing power spectra from more than one redundant baseline group, use the get_blpair_groups_from_bl_groups() method to extract certain groups.

Another useful method is fold_spectra(), which averages together \(\pm k_\parallel\) modes into a single delay spectrum as a function of \(|k_\parallel|\).

The UVPSpec class

The most relevant methods from UVPSpec are listed below. See UVPSpec Class for a full listing of all methods provided by UVPSpec.

class hera_pspec.UVPSpec

An object for storing power spectra generated by hera_pspec and a file-format for its data and meta-data.

__init__()

An object for storing power spectra and associated metadata generated by hera_pspec.

average_spectra(blpair_groups=None, time_avg=False, blpair_weights=None, error_field=None, error_weights=None, inplace=True, add_to_history='')

Average power spectra across the baseline-pair-time axis, weighted by each spectrum’s integration time.

This is an “incoherent” average, in the sense that this averages power spectra, rather than visibility data. The nsample_array and integration_array will be updated to reflect the averaging.

In the case of averaging across baseline pairs, the resultant averaged spectrum is assigned to the zeroth blpair in the group. In the case of time averaging, the time and LST arrays are assigned to the mean of the averaging window.

Note that this is designed to be separate from spherical binning in k: here we are not connecting k_perp modes to k_para modes. However, if blpairs holds groups of iso baseline separation, then this is equivalent to cylindrical binning in 3D k-space.

If you want help constructing baseline-pair groups from baseline groups, see self.get_blpair_groups_from_bl_groups.

Parameters:

blpair_groups (list) – List of list of baseline-pair group tuples or integers. All power spectra in a baseline-pair group are averaged together. If a baseline-pair exists in more than one group, a warning is raised.

Examples:

blpair_groups = [ [((1, 2), (1, 2)), ((2, 3), (2, 3))],
                  [((4, 6), (4, 6))]]

blpair_groups = [ [1002001002, 2003002003], [4006004006] ]

time_avg : bool, optional

If True, average power spectra across the time axis. Default: False.

blpair_weights : list, optional

List of float or int weights dictating the relative weight of each baseline-pair when performing the average.

This is useful for bootstrapping. This should have the same shape as blpair_groups if specified. The weights are automatically normalized within each baseline-pair group. Default: None (all baseline pairs have unity weights).

error_field: string or list, optional

If errorbars have been entered into stats_array, will do a weighted sum to shrink the error bars down to the size of the averaged data_array. Error_field strings be keys of stats_array. If list, does this for every specified key. Every stats_array key that is not specified is thrown out of the new averaged object.

error_weights: string, optional

error_weights specify which kind of errors we use for weights during averaging power spectra. The weights are defined as $w_i = 1/ sigma_i^2$, where $sigma_i$ is taken from the relevant field of stats_array.

If error_weight is set to None, which means we just use the integration time as weights. If error_weights is specified, then it also gets appended to error_field as a list. Default: None.

inplace : bool, optional

If True, edit data in self, else make a copy and return. Default: True.

add_to_history : str, optional

Added text to add to file history.

Notes

Currently, every baseline-pair in a blpair group must have the same Ntimes, unless time_avg=True. Future versions may support baseline-pair averaging of heterogeneous time arrays. This includes the scenario of repeated blpairs (e.g. in bootstrapping), which will return multiple copies of their time_array.

fold_spectra()

Average bandpowers from matching positive and negative delay bins onto a purely positive delay axis. Negative delay bins are still populated, but are filled with zero power. This is an in-place operation.

Will only work if self.folded == False, i.e. data is currently unfolded across negative and positive delay. Because this averages the data, the nsample array is multiplied by a factor of 2.

WARNING: This operation cannot be undone.

generate_noise_spectra(spw, polpair, Tsys, blpairs=None, little_h=True, form='Pk', num_steps=2000, component='real')

Generate the expected RMS noise power spectrum given a selection of spectral window, system temp. [K], and polarization. This estimate is constructed as:

P_N = scalar * (Tsys * 1e3)^2 / (integration_time) / sqrt(Nincoherent)

[mK^2 h^-3 Mpc^3]

where scalar is the cosmological and beam scalar (i.e. X2Y * Omega_eff) calculated from pspecbeam with noise_scalar = True, integration_time is in seconds and comes from self.integration_array and Nincoherent is the number of incoherent averaging samples and comes from self.nsample_array. If component is real or imag, P_N is divided by an additional factor of sqrt(2).

If the polarizations specified are pseudo Stokes pol (I, Q, U or V) then an extra factor of 2 is divided.

If form is DelSq then a factor of |k|^3 / (2pi^2) is multiplied.

If real is True, a factor of sqrt(2) is divided to account for discarding imaginary noise component.

For more details, see hera_pspec.noise.Sensitivity.calc_P_N, and Cheng et al. 2018.

The default units of P_N are mK^2 h^-3 Mpc^3 (if little_h=True and form=’Pk’), but is different if those kwargs are altered.

Parameters:

• spw (int) – Spectral window index to generate noise curve for.
• polpair (tuple or int or str) – Polarization-pair selection in form of tuple (e.g. (‘I’, ‘Q’)) or polpair int. Strings are expanded into polarization pairs, e.g. ‘XX’ becomes (‘XX,’XX’).
• Tsys (dictionary, float or array) – System temperature in Kelvin for each blpair. Key is blpair-integer, value is Tsys float or ndarray. If fed as an ndarray, shape=(Ntimes,)
• blpairs (list) – List of unique blair tuples or i12 integers to calculate noise spectrum for default is to calculate for baseline pairs.
• little_h (bool, optional) – Whether to have cosmological length units be h^-1 Mpc or Mpc Default: h^-1 Mpc
• form (str, optional) – Form of power spectra, whetehr P(k) or Delta^2(k). Valid options are ‘Pk’ or ‘DelSq’. Default: ‘Pk’.
• num_steps (int, optional) – Number of frequency bins to use in integrating power spectrum scalar in pspecbeam. Default: 2000.
• component (str) – Options=[‘real’, ‘imag’, ‘abs’]. If component is real or imag, divide by an extra factor of sqrt(2).

Returns:

P_N – Dictionary containing blpair integers as keys and float ndarrays of noise power spectra as values, with ndarrays having shape (Ntimes, Ndlys).

Return type:

dict

get_blpair_groups_from_bl_groups(blgroups, only_pairs_in_bls=False)

Get baseline pair matches from a list of baseline groups.

Parameters:

• blgroups (list) – List of baseline groups, which themselves are lists of baseline tuples or baseline i6 integers. Ex: [ [(1, 2), (2, 3), (3, 4)], [(1, 4), (2, 5)] ]
• only_pairs_in_bls (bool, optional) – If True, select only baseline-pairs whose first _and_ second baseline are both found in each baseline group. Default: False.

Returns:

blpair_groups – List of blpair groups, which themselves are lists of blpair integers.

Return type:

list

get_blpair_seps()

For each baseline-pair, get the average baseline separation in ENU frame in meters between its bl1 and bl2. Ordering matches self.get_blpairs()

Returns:blp_avg_sep – shape=(Nblpairts,) Return type:float ndarray

get_data(key, omit_flags=False)

Slice into data_array with a specified data key in the format

(spw, ((ant1, ant2), (ant3, ant4)), (pol1, pol2))

or

(spw, blpair-integer, polpair)

where spw is the spectral window integer, ant1 etc. are integers, and polpair is either a tuple of polarization strings/ints or a polarizatio-pair integer.

Parameters:

• key (tuple) – Contains the baseline-pair, spw, polpair keys
• omit_flags (bool, optional) – If True, remove time integrations (or spectra) that came from visibility data that were completely flagged across the spectral window (i.e. integration == 0).

Returns:

data – Shape (Ntimes, Ndlys)

Return type:

complex ndarray

get_dlys(key)

Get array of delays given a spectral window selection.

Parameters:key (int, or tuple with integer) – Spectral window selection Returns:dlys Return type:float ndarray, contains delays of pspectra in nanosec, given spw

get_integrations(key, omit_flags=False)

Slice into integration_array with a specified data key in the format:

(spw, ((ant1, ant2), (ant3, ant4)), (pol1, pol2))

or

(spw, blpair-integer, polpair-integer)

where spw is the spectral window integer, ant1 etc. are integers, and polpair is either a polarization pair tuple or integer.

Parameters:

• key (tuple) – Contains the baseline-pair, spw, polpair keys
• omit_flags (bool, optional) – If True, remove time integrations (or spectra) that came from visibility data that were completely flagged across the spectral window (i.e. integration == 0).

Returns:

data – Has shape (Ntimes,)

Return type:

float ndarray

get_kparas(spw, little_h=True)

Get radial (parallel) cosmological wavevectors for power spectra given an adopted cosmology and a spw selection.

Parameters:

• spw (int) – Choice of spectral window
• little_h (boolean, optional) – Whether to have cosmological length units be h^-1 Mpc or Mpc Default: h^-1 Mpc
• Returns (k_perp, k_para)
• ——-
• k_para (float ndarray) – Radial wave-vectors, shape=(Ndlys given spw)

get_kperps(spw, little_h=True)

Get transverse (perpendicular) cosmological wavevector for each baseline-pair given an adopted cosmology and a spw selection.

Parameters:

• spw (int, choice of spectral window)
• little_h (boolean, optional) – Whether to have cosmological length units be h^-1 Mpc or Mpc Default: h^-1 Mpc

Returns:

k_perp

Return type:

float ndarray, transverse wave-vectors, shape=(Nblpairs,)

get_nsamples(key, omit_flags=False)

Slice into nsample_array with a specified data key in the format

(spw, ((ant1, ant2), (ant3, ant4)), (pol1, pol2))

or

(spw, blpair-integer, polpair-integer)

where spw is the spectral window integer, ant1 etc. are integers, and polpair is either a polarization pair tuple or integer.

Parameters:

• key (tuple) – Contains the baseline-pair, spw, polpair keys
• omit_flags (bool, optional) – If True, remove time integrations (or spectra) that came from visibility data that were completely flagged across the spectral window (i.e. integration == 0).

Returns:

data

Return type:

float ndarray with shape (Ntimes,)

get_wgts(key, omit_flags=False)

Slice into wgt_array with a specified data key in the format

(spw, ((ant1, ant2), (ant3, ant4)), (pol1, pol2))

or

(spw, blpair-integer, polpair-integer)

where spw is the spectral window integer, ant1 etc. are integers, and polpair is either a polarization pair tuple or integer.

Parameters:

• key (tuple) – Contains the baseline-pair, spw, polpair keys
• omit_flags (bool, optional) – If True, remove time integrations (or spectra) that came from visibility data that were completely flagged across the spectral window (i.e. integration == 0).

Returns:

wgts – Has shape (Ntimes, Nfreqs, 2), where the last axis holds [wgt_1, wgt_2] in that order.

Return type:

float ndarray

read_hdf5(filepath, just_meta=False, spws=None, bls=None, blpairs=None, times=None, lsts=None, polpairs=None, only_pairs_in_bls=False)

Clear current UVPSpec object and load in data from an HDF5 file.

Parameters:

• filepath (str) – Path to HDF5 file to read from.

• just_meta (bool, optional) – If True, read-in metadata but ignore data, wgts and integration arrays. Default: False.

• spws (list of tuple, optional) – List of spectral window integers to select. Default: None (loads all channels).

• bls (list of i6 baseline integers or baseline tuples) – Select all baseline-pairs whose first _or_ second baseline are in the list. This changes if only_pairs_in_bls == True. Example tuple: (2, 3). Default: None (loads all bl pairs).

• blpairs (list of baseline-pair tuples or integers) – List of baseline pairs to keep. If bls is also fed, this list is concatenated onto the baseline-pair list constructed from from the bls selection.

• times (float ndarray) – Times from the time_avg_array to keep.

• lsts (float ndarray) – lsts from the lst_avg_array to keep.

• polpairs (list of polpair ints or tuples or str) – List of polarization-pair integers or tuples to keep.

  Strings are expanded into polarization pairs, and are not treated as individual polarizations, e.g. ‘XX’ becomes (‘XX,’XX’).

• only_pairs_in_bls (bool, optional) – If True, keep only baseline-pairs whose first _and_ second baseline are found in the ‘bls’ list. Default: False.

select(spws=None, bls=None, only_pairs_in_bls=True, blpairs=None, times=None, lsts=None, polpairs=None, inplace=True, run_check=True)

Select function for selecting out certain slices of the data.

Parameters:

• spws (list, optional) – List of spectral window integers to select. If None, all will be selected. Default: None.

• bls (list of i6 baseline integers or baseline tuples, optional) – Example: (2,3). Select all baseline-pairs whose first _or_ second baseline are in bls list. This changes if only_pairs_in_bls == True. If None, all baselines are selected (subject to other options). Default: None.

• only_pairs_in_bls (bool, optional) – If True, keep only baseline-pairs whose first _and_ second baseline are found in bls list. Default: True.

• blpairs (list of baseline-pair tuples or integers, optional) – List of baseline-pairs to keep. If bls is also fed, this list is concatenated onto the baseline-pair list constructed from the bls selection. If None, all valid baseline pairs are selected. Default: None.

• times (array_like, optional) – Float ndarray of times from the time_avg_array to select. If None, all times are kept. Default: None.

• lsts (array_like, optional) – Float ndarray of lsts from the lst_avg_array to select. If None, all lsts are kept. Default: None.

• polpairs (list of tuple or int or str, optional) – List of polarization-pairs to keep. If None, all polarizations are kept. Default: None.

  (Strings are expanded into polarization pairs, and are not treated as individual polarizations, e.g. ‘XX’ becomes (‘XX,’XX’).)

• inplace (bool, optional) – If True, edit and overwrite arrays in self, else make a copy of self and return. Default: True.

• run_check (bool, optional) – If True, run uvp.check() after selection

Returns:

uvp – If inplace=False, return a new UVPSpec object containing only the selected data.

Return type:

UVPSpec, optional

write_hdf5(filepath, overwrite=False, run_check=True)

Write a UVPSpec object to HDF5 file.

Parameters:

• filepath (str) – Filepath for output file.
• overwrite (bool, optional) – Whether to overwrite output file if it exists. Default: False.
• run_check (bool, optional) – Run UVPSpec validity check before writing to file. Default: True.

PSpecBeam: Primary beam models

PSpecBeam objects carry information about the primary beam, such as how the beam solid angle varies with frequency. This information is needed to rescale power spectra into cosmological units, through the computation of a ‘beam scalar’.

The main purpose of PSpecBeam objects is to provide the PSpecData class with a way of normalizing the power spectra that it produces, using the compute_pspec_scalar() method. To attach a PSpecBeam object to a PSpecData object, pass one in when you instantiate the class, e.g.

# Create PSpecBeamUV from a beamfits file
beam = hp.PSpecBeamUV('HERA_Beams.beamfits')

# Attach beam to PSpecData object
psd = hp.PSpecData(dsets=[], wgts=[], beam=beam)

Another purpose of PSpecBeam objects is to convert flux densities to temperature units using the Jy_to_mK() method, e.g.

# Apply unit conversion factor to UVData
uvd = UVData()
uvd.read_miriad(datafile) # Load data (assumed to be in Jy units)
uvd.data_array *= beam.Jy_to_mK(np.unique(uvd.freq_array))[None, None, :, None]
# (The brackets [] are needed to match the shape of uvd.data_array)

Note that PSpecBeam objects have a cosmology attached to them. If you don’t specify a cosmology (with the cosmo keyword argument), they will be instantiated with the default cosmology from hera_pspec.conversions.

There are several different types of PSpecBeam object:

PSpecBeamBase: Base class for PSpecBeam

This is the base class for all other PSpecBeam objects. It provides the generic compute_pspec_scalar() and Jy_to_mK() methods, but subclasses must provide their own power_beam_int and power_beam_sq_int methods.

class hera_pspec.pspecbeam.PSpecBeamBase(cosmo=None)

Jy_to_mK(freqs, pol='pI')

Return the multiplicative factor [mK / Jy], to convert a visibility from Jy -> mK,

factor = 1e3 * 1e-23 * c^2 / [2 * k_b * nu^2 * Omega_p(nu)]

where k_b is boltzmann constant, c is speed of light, nu is frequency and Omega_p is the integral of the unitless beam-response (steradians), and the 1e3 is the conversion from K -> mK and the 1e-23 is the conversion from Jy to cgs.

Parameters:

• freqs (float ndarray) – Contains frequencies to evaluate conversion factor [Hz].
• pol (str, optional) – Which polarization to compute the beam scalar for. ‘pI’, ‘pQ’, ‘pU’, ‘pV’, ‘XX’, ‘YY’, ‘XY’, ‘YX’ Default: ‘pI’

Returns:

factor – Contains Jy -> mK factor at each frequency.

Return type:

float ndarray

compute_pspec_scalar(lower_freq, upper_freq, num_freqs, num_steps=5000, pol='pI', taper='none', little_h=True, noise_scalar=False, exact_norm=False)

Computes the scalar function to convert a power spectrum estimate in “telescope units” to cosmological units

See arxiv:1304.4991 and HERA memo #27 for details.

Currently, only the “pI”, “XX” and “YY” polarization beams are supported. See Equations 4 and 5 of Moore et al. (2017) ApJ 836, 154 or arxiv:1502.05072 for details.

Parameters:

• lower_freq (float) – Bottom edge of frequency band over which power spectrum is being estimated. Assumed to be in Hz.
• upper_freq (float) – Top edge of frequency band over which power spectrum is being estimated. Assumed to be in Hz.
• num_freqs (int, optional) – Number of frequencies used in estimating power spectrum.
• num_steps (int, optional) – Number of steps to use when interpolating primary beams for numerical integral. Default: 5000.
• pol (str, optional) – Which polarization to compute the beam scalar for. ‘pI’, ‘pQ’, ‘pU’, ‘pV’, ‘XX’, ‘YY’, ‘XY’, ‘YX’ Default: ‘pI’
• taper (str, optional) – Whether a tapering function (e.g. Blackman-Harris) is being used in the power spectrum estimation. Default: none.
• little_h (boolean, optional) – Whether to have cosmological length units be h^-1 Mpc or Mpc. Value of h is obtained from cosmo object stored in pspecbeam. Default: h^-1 Mpc
• noise_scalar (boolean, optional) – Whether to calculate power spectrum scalar, or noise power scalar. The noise power scalar only differs in that the Bpp_over_BpSq term just because 1_over_Bp. See Pober et al. 2014, ApJ 782, 66.
• exact_norm (boolean, optional) – returns only X2Y for scalar if True, else uses the existing framework involving antenna beam and spectral tapering factors. Default: False.

Returns:

scalar – [int dnu (Omega_PP / Omega_P^2) ( B_PP / B_P^2 ) / (X^2 Y)]^-1 Units: h^-3 Mpc^3 or Mpc^3.

Return type:

float

get_Omegas(polpairs)

Get OmegaP and OmegaPP across beam_freqs for requested polarization pairs.

Parameters:polpairs (list) – List of polarization-pair tuples or integers. Returns:

• OmegaP (array_like) – Array containing power_beam_int, shape: (Nbeam_freqs, Npols).
• OmegaPP (array_like) – Array containing power_sq_beam_int, shape: (Nbeam_freqs, Npols).

PSpecBeamUV: Beams from a UVBeam object

This class allows you to use any beam that is supported by the UVBeam class in the pyuvdata package. These usually contain Healpix-pixelated beams as a function of frequency and polarization.

To create a beam that uses this format, simply create a new PSpecBeamUV instance with the name of a beamfits file that is supported by UVBeam, e.g.

beam = hp.PSpecBeamUV('HERA_Beams.beamfits')

class hera_pspec.pspecbeam.PSpecBeamUV(uvbeam, cosmo=None)

beam_normalized_response(pol='pI', freq=None, x_orientation=None)

Outputs beam response for given polarization as a function of pixels on the sky and input frequencies. The response needs to be peak normalized, and is read in from Healpix coordinates. Uses interp_freq function from uvbeam for interpolation of beam response over given frequency values.

Parameters:

• pol (str, optional) – Which polarization to compute the beam response for. ‘pI’, ‘pQ’, ‘pU’, ‘pV’, ‘XX’, ‘YY’, ‘XY’, ‘YX’ The output shape is (Nfreq, Npixels) Default: ‘pI’
• freq (array, optional) – Frequencies [Hz] to interpolate onto.
• x_orientation (str, optional) – Orientation in cardinal direction east or north of X dipole. Default keeps polarization in X and Y basis.

Returns:

• beam_res (float, array-like) – Beam response as a function healpix indices and frequency.
• omega (float, array-like) – Beam solid angle as a function of frequency
• nside (int, scalar) – used to compute resolution

power_beam_int(pol='pI')

Computes the integral of the beam over solid angle to give a beam area (in str) as a function of frequency. Uses function in pyuvdata.

See Equations 4 and 5 of Moore et al. (2017) ApJ 836, 154 or arxiv:1502.05072 for details.

Parameters:pol (str, optional) – Which polarization to compute the beam scalar for. ‘pI’, ‘pQ’, ‘pU’, ‘pV’, ‘XX’, ‘YY’, ‘XY’, ‘YX’ Default: ‘pI’ Returns:primary_beam_area – Scalar integral over beam solid angle. Return type:float, array-like

power_beam_sq_int(pol='pI')

Computes the integral of the beam**2 over solid angle to give a beam**2 area (in str) as a function of frequency. Uses function in pyuvdata.

See Equations 4 and 5 of Moore et al. (2017) ApJ 836, 154 or arxiv:1502.05072 for details.

Parameters:pol (str, optional) – Which polarization to compute the beam scalar for. ‘pI’, ‘pQ’, ‘pU’, ‘pV’, ‘XX’, ‘YY’, ‘XY’, ‘YX’ Default: ‘pI’ Returns:primary_beam_area Return type:float, array-like

Jy_to_mK(freqs, pol='pI')

Return the multiplicative factor [mK / Jy], to convert a visibility from Jy -> mK,

factor = 1e3 * 1e-23 * c^2 / [2 * k_b * nu^2 * Omega_p(nu)]

where k_b is boltzmann constant, c is speed of light, nu is frequency and Omega_p is the integral of the unitless beam-response (steradians), and the 1e3 is the conversion from K -> mK and the 1e-23 is the conversion from Jy to cgs.

Parameters:

• freqs (float ndarray) – Contains frequencies to evaluate conversion factor [Hz].
• pol (str, optional) – Which polarization to compute the beam scalar for. ‘pI’, ‘pQ’, ‘pU’, ‘pV’, ‘XX’, ‘YY’, ‘XY’, ‘YX’ Default: ‘pI’

Returns:

factor – Contains Jy -> mK factor at each frequency.

Return type:

float ndarray

compute_pspec_scalar(lower_freq, upper_freq, num_freqs, num_steps=5000, pol='pI', taper='none', little_h=True, noise_scalar=False, exact_norm=False)

Computes the scalar function to convert a power spectrum estimate in “telescope units” to cosmological units

See arxiv:1304.4991 and HERA memo #27 for details.

Currently, only the “pI”, “XX” and “YY” polarization beams are supported. See Equations 4 and 5 of Moore et al. (2017) ApJ 836, 154 or arxiv:1502.05072 for details.

Parameters:

• lower_freq (float) – Bottom edge of frequency band over which power spectrum is being estimated. Assumed to be in Hz.
• upper_freq (float) – Top edge of frequency band over which power spectrum is being estimated. Assumed to be in Hz.
• num_freqs (int, optional) – Number of frequencies used in estimating power spectrum.
• num_steps (int, optional) – Number of steps to use when interpolating primary beams for numerical integral. Default: 5000.
• pol (str, optional) – Which polarization to compute the beam scalar for. ‘pI’, ‘pQ’, ‘pU’, ‘pV’, ‘XX’, ‘YY’, ‘XY’, ‘YX’ Default: ‘pI’
• taper (str, optional) – Whether a tapering function (e.g. Blackman-Harris) is being used in the power spectrum estimation. Default: none.
• little_h (boolean, optional) – Whether to have cosmological length units be h^-1 Mpc or Mpc. Value of h is obtained from cosmo object stored in pspecbeam. Default: h^-1 Mpc
• noise_scalar (boolean, optional) – Whether to calculate power spectrum scalar, or noise power scalar. The noise power scalar only differs in that the Bpp_over_BpSq term just because 1_over_Bp. See Pober et al. 2014, ApJ 782, 66.
• exact_norm (boolean, optional) – returns only X2Y for scalar if True, else uses the existing framework involving antenna beam and spectral tapering factors. Default: False.

Returns:

scalar – [int dnu (Omega_PP / Omega_P^2) ( B_PP / B_P^2 ) / (X^2 Y)]^-1 Units: h^-3 Mpc^3 or Mpc^3.

Return type:

float

get_Omegas(polpairs)

Get OmegaP and OmegaPP across beam_freqs for requested polarization pairs.

Parameters:polpairs (list) – List of polarization-pair tuples or integers. Returns:

• OmegaP (array_like) – Array containing power_beam_int, shape: (Nbeam_freqs, Npols).
• OmegaPP (array_like) – Array containing power_sq_beam_int, shape: (Nbeam_freqs, Npols).

PSpecBeamGauss: Simple Gaussian beam model

A Gaussian beam type is provided for simple testing purposes. You can specify a beam FWHM that is constant in frequency, for the I (pseudo-Stokes I) polarization channel only.

For example, to specify a Gaussian beam with a constant FWHM of 0.1 radians, defined over a frequency interval of [100, 200] MHz:

# Each beam is defined over a frequency interval:
beam_freqs = np.linspace(100e6, 200e6, 200) # in Hz

# Create a new Gaussian beam object with full-width at half-max. of 0.1 radians
beam_gauss = hp.PSpecBeamGauss(fwhm=0.1, beam_freqs=beam_freqs)

class hera_pspec.pspecbeam.PSpecBeamGauss(fwhm, beam_freqs, cosmo=None)

__init__(fwhm, beam_freqs, cosmo=None)

Object to store a simple (frequency independent) Gaussian beam in a PspecBeamBase object.

Parameters:

• fwhm (float) – Full width half max of the beam, in radians.
• beam_freqs (float, array-like) – Frequencies over which this Gaussian beam is to be created. Units assumed to be Hz.
• cosmo (conversions.Cosmo_Conversions object, optional) – Cosmology object. Uses the default cosmology object if not specified. Default: None.

power_beam_int(pol='pI')

Computes the integral of the beam over solid angle to give a beam area (in sr). Uses analytic formula that the answer is 2 * pi * fwhm**2 / 8 ln 2.

Trivially this returns an array (i.e., a function of frequency), but the results are frequency independent.

See Equations 4 and 5 of Moore et al. (2017) ApJ 836, 154 or arxiv:1502.05072 for details.

Parameters:pol (str, optional) – Which polarization to compute the beam scalar for. ‘pI’, ‘pQ’, ‘pU’, ‘pV’, ‘XX’, ‘YY’, ‘XY’, ‘YX’ Default: ‘pI’ Returns:primary_beam_area – Primary beam area. Return type:float, array-like

power_beam_sq_int(pol='pI')

Computes the integral of the beam**2 over solid angle to give a beam area (in str). Uses analytic formula that the answer is pi * fwhm**2 / 8 ln 2.

Trivially this returns an array (i.e., a function of frequency), but the results are frequency independent.

See Equations 4 and 5 of Moore et al. (2017) ApJ 836, 154 or arxiv:1502.05072 for details.

Parameters:pol (str, optional) – Which polarization to compute the beam scalar for. ‘pI’, ‘pQ’, ‘pU’, ‘pV’, ‘XX’, ‘YY’, ‘XY’, ‘YX’ Default: ‘pI’ Returns:primary_beam_area – Primary beam area. Return type:float, array-like

PSpecData Class

PSpecData takes UVData objects and calculates delay power spectra from them.

class hera_pspec.PSpecData(dsets=[], wgts=None, dsets_std=None, labels=None, beam=None, cals=None, cal_flag=True)

C_model(key, model='empirical', time_index=None, known_cov=None, include_extension=False)

Return a covariance model having specified a key and model type. Note: Time-dependent flags that differ from frequency channel-to-channel can create spurious spectral structure. Consider factorizing the flags with self.broadcast_dset_flags() before using model=’empirical’.

Parameters:

• key (tuple) – Tuple containing indices of dataset and baselines. The first item specifies the index (ID) of a dataset in the collection, while subsequent indices specify the baseline index, in _key2inds format.
• model (string, optional) – Type of covariance model to calculate, if not cached. Options=[‘empirical’, ‘dsets’, ‘autos’, (other model names in known_cov)] How the covariances of the input data should be estimated. In ‘dsets’ mode, error bars are estimated from user-provided per baseline and per channel standard deivations. If ‘empirical’ is set, then error bars are estimated from the data by averaging the channel-channel covariance of each baseline over time and then applying the appropriate linear transformations to these frequency-domain covariances. If ‘autos’ is set, the covariances of the input data over a baseline is estimated from the autocorrelations of the two antennas over channel bandwidth and integration time.
• time_index (integer, compute covariance at specific time-step in dset) – supported if mode == ‘dsets’ or ‘autos’
• known_cov (dicts of covariance matrices) – Covariance matrices that are imported from a outer dict instead of using data stored or calculated inside the PSpecData object. known_cov could be initialized when using PSpecData.pspec() method. See PSpecData.pspec() for more details.
• include_extension (bool (optional)) – default=False If True, extend spw to include filtering window extensions.

Returns:

C – Covariance model for the specified key.

Return type:

ndarray, (spw_Nfreqs, spw_Nfreqs)

I(key)

Return identity covariance matrix.

Parameters:key (tuple) – Tuple containing indices of dataset and baselines. The first item specifies the index (ID) of a dataset in the collection, while subsequent indices specify the baseline index, in _key2inds format. Returns:I – Identity covariance matrix, dimension (Nfreqs, Nfreqs). Return type:array_like

Jy_to_mK(beam=None)

Convert internal datasets from a Jy-scale to mK scale using a primary beam model if available. Note that if you intend to rephase_to_dset(), Jy to mK conversion must be done after that step.

Parameters:beam (PSpecBeam object) – Beam object.

R(key)

Return the data-weighting matrix R, which is a product of data covariance matrix (I or C^-1), diagonal flag matrix (Y) and diagonal tapering matrix (T):

R = sqrt(T^t) sqrt(Y^t) K sqrt(Y) sqrt(T)

where T is a diagonal matrix holding the taper and Y is a diagonal matrix holding flag weights. The K matrix comes from either I or iC or a dayenu depending on self.data_weighting, T is informed by self.taper and Y is taken from self.Y().

Right now, the data covariance can be identity (‘I’), C^-1 (‘iC’), or dayenu weighting ‘dayenu’.

Parameters:key (tuple) – Tuple containing indices of dataset and baselines. The first item specifies the index (ID) of a dataset in the collection, while subsequent indices specify the baseline index, in _key2inds format.

Y(key)

Return the weighting (diagonal) matrix, Y. This matrix is calculated by taking the logical AND of flags across all times given the dset-baseline-pol specification in ‘key’, converted into a float, and inserted along the diagonal of an spw_Nfreqs x spw_Nfreqs matrix.

The logical AND step implies that all time-dependent flagging patterns are automatically broadcasted across all times. This broadcasting follows the principle that, for each freq channel, if at least a single time is unflagged, then the channel is treated as unflagged for all times. Power spectra from certain times, however, can be given zero weight by setting the nsample array to be zero at those times (see self.broadcast_dset_flags).

Parameters:key (tuple) – Tuple containing indices of dataset and baselines. The first item specifies the index (ID) of a dataset in the collection, while subsequent indices specify the baseline index, in _key2inds format. Returns:Y – spw_Nfreqs x spw_Nfreqs diagonal matrix holding AND of flags across all times for each freq channel. Return type:array_like

add(dsets, wgts, labels=None, dsets_std=None, cals=None, cal_flag=True)

Add a dataset to the collection in this PSpecData object.

Parameters:

• dsets (UVData or list or dict) – UVData object or list of UVData objects containing data to add to the collection.
• wgts (UVData or list or dict) – UVData object or list of UVData objects containing weights to add to the collection. Must be the same length as dsets. If a weight is set to None, the flags of the corresponding dset are used.
• labels (list of str) – An ordered list of names/labels for each dataset, if dsets was specified as a list. If dsets was specified as a dict, the keys of that dict will be used instead.
• dsets_std (UVData or list or dict) – Optional UVData object or list of UVData objects containing the standard deviations (real and imaginary) of data to add to the collection. If dsets is a dict, will assume dsets_std is a dict and if dsets is a list, will assume dsets_std is a list.
• cals (UVCal or list, optional) – UVCal objects to apply to data.
• cal_flag (bool, optional) – If True, propagate flags from calibration into data

broadcast_dset_flags(spw_ranges=None, time_thresh=0.2, unflag=False)

For each dataset in self.dset, update the flag_array such that the flagging patterns are time-independent for each baseline given a selection for spectral windows.

For each frequency pixel in a selected spw, if the fraction of flagged times exceeds time_thresh, then all times are flagged. If it does not, the specific integrations which hold flags in the spw are flagged across all frequencies in the spw.

Additionally, one can also unflag the flag_array entirely if desired.

Note: although technically allowed, this function may give unexpected results if multiple spectral windows in spw_ranges have frequency overlap.

Note: it is generally not recommended to set time_thresh > 0.5, which could lead to substantial amounts of data being flagged.

Parameters:

• spw_ranges (list of tuples) – list of len-2 spectral window tuples, specifying the start (inclusive) and stop (exclusive) index of the frequency array for each spw. Default is to use the whole band.
• time_thresh (float) – Fractional threshold of flagged pixels across time needed to flag all times per freq channel. It is not recommend to set this greater than 0.5.
• unflag (bool) – If True, unflag all data in the spectral window.

check_key_in_dset(key, dset_ind)

Check ‘key’ exists in the UVData object self.dsets[dset_ind]

Parameters:

• key (tuple) – if length 1: assumed to be polarization number or string elif length 2: assumed to be antenna-number tuple (ant1, ant2) elif length 3: assumed ot be antenna-number-polarization tuple (ant1, ant2, pol)
• dset_ind (int, the index of the dataset to-be-checked)

Returns:

exists – True if the key exists, False otherwise

Return type:

bool

clear_cache(keys=None)

Clear stored matrix data (or some subset of it).

Parameters:keys (list of tuples, optional) – List of keys to remove from matrix cache. If None, all keys will be removed. Default: None.

cov_p_hat(M, q_cov)

Covariance estimate between two different band powers p_alpha and p_beta given by M_{alpha i} M^*_{beta,j} C_q^{ij} where C_q^{ij} is the q-covariance.

Parameters:

• M (array_like) – Normalization matrix, M.
• q_cov (array_like) – covariance between bandpowers in q_alpha and q_beta

cov_q_hat(key1, key2, model='empirical', exact_norm=False, pol=False, time_indices=None)

Compute the un-normalized covariance matrix for q_hat for a given pair of visibility vectors. Returns the following matrix:

Cov(hat{q}_a,hat{q}_b)

!!!Only supports covariance between same power-spectrum estimates!!! (covariance between pair of baselines with the same pair of baselines) !!!Assumes that both baselines used in power-spectrum estimate !!!have independent noise relizations!!!

#updated to be a multi-time wrapper to get_unnormed_V

Parameters:

• key1, key2 (tuples or lists of tuples) – Tuples containing indices of dataset and baselines for the two input datavectors. If a list of tuples is provided, the baselines in the list will be combined with inverse noise weights.
• exact_norm (boolean, optional) – Exact normalization (see HERA memo #44, Eq. 11 and documentation of q_hat for details).
• pol (str/int/bool, optional) – Polarization parameter to be used for extracting the correct beam. Used only if exact_norm is True.
• model (string, optional) – Type of covariance model to calculate, if not cached. Options=[‘empirical’, ‘dsets’, ‘autos’] How the covariances of the input data should be estimated. In ‘dsets’ mode, error bars are estimated from user-provided per baseline and per channel standard deivations. If ‘empirical’ is set, then error bars are estimated from the data by averaging the channel-channel covariance of each baseline over time and then applying the appropriate linear transformations to these frequency-domain covariances. If ‘autos’ is set, the covariances of the input data over a baseline is estimated from the autocorrelations of the two antennas over channel bandwidth and integration time.
• time_indices (list of indices of times to include or just a single time.)
• default is None -> compute covariance for all times.

Returns:

cov_q_hat – Matrix with covariances between un-normalized band powers (Ntimes, Nfreqs, Nfreqs)

Return type:

array_like

cross_covar_model(key1, key2, model='empirical', time_index=None, conj_1=False, conj_2=True, known_cov=None, include_extension=False)

Return a covariance model having specified a key and model type. Note: Time-dependent flags that differ from frequency channel-to-channel can create spurious spectral structure. Consider factorizing the flags with self.broadcast_dset_flags() before using model=’time_average’.

Parameters:

• key1, key2 (tuples) – Tuples containing indices of dataset and baselines. The first item specifies the index (ID) of a dataset in the collection, while subsequent indices specify the baseline index, in _key2inds format.
• model (string, optional) – Type of covariance model to calculate, if not cached. Options=[‘empirical’, ‘dsets’, ‘autos’, (other model names in known_cov)] How the covariances of the input data should be estimated. In ‘dsets’ mode, error bars are estimated from user-provided per baseline and per channel standard deivations. If ‘empirical’ is set, then error bars are estimated from the data by averaging the channel-channel covariance of each baseline over time and then applying the appropriate linear transformations to these frequency-domain covariances. If ‘autos’ is set, the covariances of the input data over a baseline is estimated from the autocorrelations of the two antennas over channel bandwidth and integration time.
• time_index (integer, compute covariance at specific time-step)
• conj_1 (boolean, optional) – Whether to conjugate first copy of data in covar or not. Default: False
• conj_2 (boolean, optional) – Whether to conjugate second copy of data in covar or not. Default: True
• known_cov (dicts of covariance matrices) – Covariance matrices that are imported from a outer dict instead of using data stored or calculated inside the PSpecData object. known_cov could be initialized when using PSpecData.pspec() method. See PSpecData.pspec() for more details.
• include_extension (bool (optional)) – default=False If True, extend spw to include filtering window extensions.

Returns:

cross_covar – Cross covariance model for the specified key.

Return type:

ndarray, (spw_Nfreqs, spw_Nfreqs)

delays()

Return an array of delays, tau, corresponding to the bins of the delay power spectrum output by pspec() using self.spw_range to specify the spectral window.

Returns:delays – Delays, tau. Units: ns. Return type:array_like

dset_idx(dset)

Return the index of a dataset, regardless of whether it was specified as an integer of a string.

Parameters:dset (int or str) – Index or name of a dataset belonging to this PSpecData object. Returns:dset_idx – Index of dataset. Return type:int

dx(key, include_extension=False)

Get standard deviation of data for given dataset and baseline as pecified in standard key format.

Parameters:

• key (tuple) – Tuple containing datset ID and baseline index. The first element of the tuple is the dataset index (or label), and the subsequent elements are the baseline ID.
• include_extension (bool (optional)) – default=False If True, extend spw to include filtering window extensions.

Returns:

dx – Array of std data from the requested UVData dataset and baseline.

Return type:

array_like

get_G(key1, key2, exact_norm=False, pol=False)

Calculates

G_ab = (1/2) Tr[R_1 Q^alt_a R_2 Q^alt_b],

which is needed for normalizing the power spectrum (see HERA memo #44).

Note that in the limit that R_1 = R_2 = C^-1, this reduces to the Fisher matrix

F_ab = 1/2 Tr [C^-1 Q^alt_a C^-1 Q^alt_b] (arXiv:1502.06016, Eq. 17)

Parameters:

• key1, key2 (tuples or lists of tuples) – Tuples containing indices of dataset and baselines for the two input datavectors. If a list of tuples is provided, the baselines in the list will be combined with inverse noise weights.
• exact_norm (boolean, optional) – Exact normalization (see HERA memo #44, Eq. 11 and documentation of q_hat for details).
• pol (str/int/bool, optional) – Polarization parameter to be used for extracting the correct beam. Used only if exact_norm is True.

Returns:

G – Fisher matrix, with dimensions (Nfreqs, Nfreqs).

Return type:

array_like, complex

get_H(key1, key2, sampling=False, exact_norm=False, pol=False)

Calculates the response matrix H of the unnormalized band powers q to the true band powers p, i.e.,

<q_a> = sum_b H_{ab} p_b

This is given by

H_ab = (1/2) Tr[R_1 Q_a^alt R_2 Q_b]

(See HERA memo #44). As currently implemented, this approximates the primary beam as frequency independent. Under this approximation, the our H_ab is defined using the equation above except we have Q^tapered rather than Q_b, where

overline{Q}^{tapered,beta} = e^{i 2pi eta_beta (nu_i - nu_j)} gamma(nu_i) gamma(nu_j)

where gamma is the tapering function. Again, see HERA memo #44 for details.

The sampling option determines whether one is assuming that the output points are integrals over k bins or samples at specific k values. The effect is to add a dampening of widely separated frequency correlations with the addition of the term

sinc(pi Delta eta (nu_i - nu_j))

Note that numpy uses the engineering definition of sinc, where sinc(x) = sin(pi x) / (pi x), whereas in the line above and in all of our documentation, we use the physicist definition where sinc(x) = sin(x) / x

Note that in the limit that R_1 = R_2 = C^-1 and Q_a is used instead of Q_a^alt, this reduces to the Fisher matrix

F_ab = 1/2 Tr [C^-1 Q_a C^-1 Q_b] (arXiv:1502.06016, Eq. 17)

This function uses the state of self.taper in constructing H. See PSpecData.pspec for details.

Parameters:

• key1, key2 (tuples or lists of tuples) – Tuples containing indices of dataset and baselines for the two input datavectors. If a list of tuples is provided, the baselines in the list will be combined with inverse noise weights.
• sampling (boolean, optional) – Whether to sample the power spectrum or to assume integrated bands over wide delay bins. Default: False
• exact_norm (boolean, optional) – Exact normalization (see HERA memo #44, Eq. 11 and documentation of q_hat for details).
• pol (str/int/bool, optional) – Polarization parameter to be used for extracting the correct beam. Used only if exact_norm is True.

Returns:

H – Dimensions (Nfreqs, Nfreqs).

Return type:

array_like, complex

get_MW(G, H, mode='I', band_covar=None, exact_norm=False, rcond=1e-15)

Construct the normalization matrix M and window function matrix W for the power spectrum estimator. These are defined through Eqs. 14-16 of arXiv:1502.06016:

hat{p} = M hat{q} <hat{p}> = W p W = M H,

where p is the true band power and H is the response matrix (defined above in get_H) of unnormalized bandpowers to normed bandpowers.

Several choices for M are supported:

‘I’: Set M to be diagonal (e.g. HERA Memo #44) ‘H^-1’: Set M = H^-1, the (pseudo)inverse response matrix. ‘V^-1/2’: Set M = V^-1/2, the root-inverse response matrix (using SVD).

These choices will be supported very soon:

‘L^-1’: Set M = L^-1, Cholesky decomposition.

As written, the window functions will not be correclty normalized; it needs to be adjusted by the pspec scalar for it to be approximately correctly normalized. If the beam is being provided, this will be done in the pspec function.

Parameters:

• G (array_like) – Denominator matrix for the bandpowers, with dimensions (Nfreqs, Nfreqs).
• H (array_like) – Response matrix for the bandpowers, with dimensions (Nfreqs, Nfreqs).
• mode (str, optional) – Definition to use for M. Must be one of the options listed above. Default: ‘I’.
• band_covar (array_like, optional) – Covariance matrix of the unnormalized bandpowers (i.e., q). Used only if requesting the V^-1/2 normalization. Use get_unnormed_V to get the covariance to put in here, or provide your own array. Default: None
• exact_norm (boolean, optional) – Exact normalization (see HERA memo #44, Eq. 11 and documentation of q_hat for details). Currently, this is supported only for mode I
• rcond (float, optional) – rcond parameter of np.linalg.pinv for truncating near-zero eigenvalues

Returns:

• M (array_like) – Normalization matrix, M. (If G was passed in as a dict, a dict of array_like will be returned.)
• W (array_like) – Window function matrix, W. (If G was passed in as a dict, a dict of array_like will be returned.)

get_Q(mode)

Computes Q_alt(i,j), which is the exponential part of the response of the data covariance to the bandpower (dC/dP_alpha).

Note: This function is not being used right now, since get_q_alt and get_Q are essentially the same functions. However, since we want to attempt non-uniform bins, we do intend to use get_Q (which uses physical units, and hence there is not contraint of uniformly spaced data).

Parameters:mode (int) – Central wavenumber (index) of the bandpower, p_alpha. Returns:Q_alt – Exponential part of Q (HERA memo #44, Eq. 11). Return type:array_like

get_Q_alt(mode, allow_fft=True, include_extension=False)

Response of the covariance to a given bandpower, dC / dp_alpha, EXCEPT without the primary beam factors. This is Q_alt as defined in HERA memo #44, so it’s not dC / dp_alpha, strictly, but is just the part that does the Fourier transforms.

Assumes that Q will operate on a visibility vector in frequency space. In the limit that self.spw_Ndlys equals self.spw_Nfreqs, this will produce a matrix Q that performs a two-sided FFT and extracts a particular Fourier mode.

(Computing x^t Q y is equivalent to Fourier transforming x and y separately, extracting one element of the Fourier transformed vectors, and then multiplying them.)

When self.spw_Ndlys < self.spw_Nfreqs, the effect is similar except the delay bins need not be in the locations usually mandated by the FFT algorithm.

Parameters:

• mode (int) – Central wavenumber (index) of the bandpower, p_alpha.
• allow_fft (boolean, optional) – If set to True, allows a shortcut FFT method when the number of delay bins equals the number of delay channels. Default: True
• include_extension (If True, return a matrix that is spw_Nfreq x spw_Nfreq)
• (required if using partial C_{ij} / partial p_lpha since C_{ij} is
• (spw_Nfreq x spw_Nfreq).

Returns:

Q – Response matrix for bandpower p_alpha.

Return type:

array_like

get_analytic_covariance(key1, key2, M=None, exact_norm=False, pol=False, model='empirical', known_cov=None)

Calculates the auto-covariance matrix for both the real and imaginary parts of bandpowers (i.e., the q vectors and the p vectors).

Define:

Real part of q_a = (1/2) (q_a + q_a^*) Imaginary part of q_a = (1/2i) (q_a - q_a^dagger) Real part of p_a = (1/2) (p_a + p_a^dagger) Imaginary part of p_a = (1/2i) (p_a - p_a^dagger)

\[E^{12,a} = (1/2) R_1 Q^a R_2 C^{12} = <x1 x2^\dagger> - <x1><x2^\dagger> P^{12} = <x1 x2> - <x1><x2> S^{12} = <x1^* x2^*> - <x1^*> <x2^*> p_a = M_{ab} q_b\]

Then:

The variance of (1/2) (q_a + q_a^dagger):

\[(1/4){ (<q_a q_a> - <q_a><q_a>) + 2(<q_a q_a^\dagger> - <q_a><q_a^\dagger>) + (<q_a^\dagger q_a^\dagger> - <q_a^\dagger><q_a^\dagger>) }\]

The variance of (1/2i) (q_a - q_a^dagger):

\[(-1/4){ (<q_a q_a> - <q_a><q_a>) - 2(<q_a q_a^\dagger> - <q_a><q_a^\dagger>) + (<q_a^\dagger q_a^\dagger> - <q_a^\dagger><q_a^\dagger>) }\]

The variance of (1/2) (p_a + p_a^dagger):

\[(1/4) { M_{ab} M_{ac} (<q_b q_c> - <q_b><q_c>) + M_{ab} M_{ac}^* (<q_b q_c^\dagger> - <q_b><q_c^\dagger>) + M_{ab}^* M_{ac} (<q_b^\dagger q_c> - <q_b^\dagger><q_c>) + M_{ab}^* M_{ac}^* (<q_b^\dagger q_c^\dagger> - <q_b^\dagger><q_c^\dagger>) }\]

The variance of (1/2i) (p_a - p_a^dagger):

\[(-1/4) { M_{ab} M_{ac} (<q_b q_c> - <q_b><q_c>) - M_{ab} M_{ac}^* (<q_b q_c^\dagger> - <q_b><q_c^\dagger>) - M_{ab}^* M_{ac} (<q_b^\dagger q_c> - <q_b^\dagger><q_c>) + M_{ab}^* M_{ac}^* (<q_b^\dagger q_c^\dagger> - <q_b^\dagger><q_c^\dagger>) }\]

where

\[<q_a q_b> - <q_a><q_b> = tr(E^{12,a} C^{21} E^{12,b} C^{21}) + tr(E^{12,a} P^{22} E^{21,b*} S^{11}) <q_a q_b^\dagger> - <q_a><q_b^\dagger> = tr(E^{12,a} C^{22} E^{21,b} C^{11}) + tr(E^{12,a} P^{21} E^{12,b *} S^{21}) <q_a^\dagger q_b^\dagger> - <q_a^\dagger><q_b^\dagger> = tr(E^{21,a} C^{12} E^{21,b} C^{12}) + tr(E^{21,a} P^{11} E^{12,b *} S^{22})\]

Note that

\[E^{12,a}_{ij}.conj = E^{21,a}_{ji}\]

This function estimates C^1, C^2, P^{12}, and S^{12} empirically by default. (So while the pointy brackets <…> should in principle be ensemble averages, in practice the code performs averages in time.)

Note: Time-dependent flags that differ from frequency channel-to-channel can create spurious spectral structure. Consider factorizing the flags with self.broadcast_dset_flags() before using model=’time_average’

Parameters:

• key1, key2 (tuples or lists of tuples) – Tuples containing indices of dataset and baselines for the two input datavectors. If a list of tuples is provided, the baselines in the list will be combined with inverse noise weights.

• M (array_like) – Normalization matrix, M. Ntimes x Ndlys x Ndlys

• exact_norm (boolean) – If True, beam and spectral window factors are taken in the computation of Q_matrix (dC/dp = Q, and not Q_alt) (HERA memo #44, Eq. 11). Q matrix, for each delay mode, is weighted by the integral of beam over theta,phi. Therefore the output power spectra is, by construction, normalized. If True, it returns normalized power spectrum, except for X2Y term. If False, Q_alt is used (HERA memo #44, Eq. 16), and the power spectrum is normalized separately.

• pol (str/int/bool, optional) – Polarization parameter to be used for extracting the correct beam. Used only if exact_norm is True.

• model (string, optional) – Type of covariance model to use. if not cached. Options=[‘empirical’, ‘dsets’, ‘autos’, ‘foreground_dependent’, (other model names in known_cov)].

  In dsets mode, error bars are estimated from user-provided per baseline and per channel standard deivations.

  In empirical mode, error bars are estimated from the data by averaging the channel-channel covariance of each baseline over time and then applying the appropriate linear transformations to these frequency-domain covariances.

  In autos mode, the covariances of the input data over a baseline is estimated from the autocorrelations of the two antennas forming the baseline across channel bandwidth and integration time.

  In foreground_dependent mode, it involves using auto-correlation amplitudes to model the input noise covariance and visibility outer products to model the input systematics covariance.

  When model is chosen as autos or dsets, only C^{11} and C^{22} are accepted as non-zero values, and the two matrices are also expected to be diagonal, thus only <q_a q_b^dagger> - <q_a><q_b^dagger> = tr[ E^{12,a} C^{22} E^{21,b} C^{11} ] exists in the covariance terms of q vectors.

  When model is chosen as foreground_dependent, we further include the signal-noise coupling term besides the noise in the output covariance. Still only <q_a q_b^dagger> - <q_a><q_b^dagger> is non-zero, while it takes a form of tr[ E^{12,a} Cn^{22} E^{21,b} Cn^{11} + E^{12,a} Cs^{22} E^{21,b} Cn^{11} + E^{12,a} Cn^{22} E^{21,b} Cs^{11} ], where Cn is just Cautos, the input noise covariance estimated by the auto-correlation amplitudes (by calling C_model(model=’autos’)), and Cs uses the outer product of input visibilities to model the covariance on systematics.

  To construct a symmetric and unbiased covariance matrix, we choose Cs^{11}_{ij} = Cs^{22}_{ij} = 1/2 * [ x1_i x2_j^{*} + x2_i x1_j^{*} ], which preserves the property Cs_{ij}^* = Cs_{ji}.

• known_cov (dicts of covariance matrices) – Covariance matrices that are not calculated internally from data.

Returns:

V – Bandpower covariance, with dimension (Ntimes, spw_Ndlys, spw_Ndlys).

Return type:

array_like, complex

get_integral_beam(pol=False)

Computes the integral containing the spectral beam and tapering function in Q_alpha(i,j).

Parameters:pol (str/int/bool, optional) – Which beam polarization to use. If the specified polarization doesn’t exist, a uniform isotropic beam (with integral 4pi for all frequencies) is assumed. Default: False (uniform beam). Returns:integral_beam – integral containing the spectral beam and tapering. Return type:array_like

get_spw(include_extension=False)

Get self.spw_range with or without spw extension (self.filter_extension)

Parameters:include_extension (bool) – If True, include self.filter_extension in spw_range Returns:spectral_window – In (start_chan, end_chan). Return type:tuple

get_unnormed_E(key1, key2, exact_norm=False, pol=False)

Calculates a series of unnormalized E matrices, such that

q_a = x_1^* E^{12,a} x_2

so that

E^{12,a} = (1/2) R_1 Q^a R_2.

In principle, this could be used to actually estimate q-hat. In other words, we could call this function to get E, and then sandwich it with two data vectors to get q_a. However, this should be slower than the implementation in q_hat. So this should only be used as a helper method for methods such as get_unnormed_V.

Note for the future: There may be advantages to doing the matrix multiplications separately

Parameters:

• key1, key2 (tuples or lists of tuples) – Tuples containing indices of dataset and baselines for the two input datavectors. If a list of tuples is provided, the baselines in the list will be combined with inverse noise weights.
• exact_norm (boolean, optional) – Exact normalization (see HERA memo #44, Eq. 11 and documentation of q_hat for details).
• pol (str/int/bool, optional) – Polarization parameter to be used for extracting the correct beam. Used only if exact_norm is True.

Returns:

E – Set of E matrices, with dimensions (Ndlys, Nfreqs, Nfreqs).

Return type:

array_like, complex

get_unnormed_V(key1, key2, model='empirical', exact_norm=False, pol=False, time_index=None)

Calculates the covariance matrix for unnormed bandpowers (i.e., the q vectors). If the data were real and x_1 = x_2, the expression would be

\[V_ab = 2 tr(C E_a C E_b), where E_a = (1/2) R Q^a R\]

When the data are complex, the expression becomes considerably more complicated. Define

\[E^{12,a} = (1/2) R_1 Q^a R_2 C^1 = <x1 x1^\dagger> - <x1><x1^\dagger> C^2 = <x2 x2^\dagger> - <x2><x2^\dagger> P^{12} = <x1 x2> - <x1><x2> S^{12} = <x1^* x2^*> - <x1^*> <x2^*>\]

Then

\[V_ab = tr(E^{12,a} C^2 E^{21,b} C^1) + tr(E^{12,a} P^{21} E^{12,b *} S^{21})\]

Note that

\[E^{12,a}_{ij}.conj = E^{21,a}_{ji}\]

This function estimates C^1, C^2, P^{12}, and S^{12} empirically by default. (So while the pointy brackets <…> should in principle be ensemble averages, in practice the code performs averages in time.)

Empirical covariance estimates are in principle a little risky, as they can potentially induce signal loss. This is probably ok if we are just looking intending to look at V. It is most dangerous when C_emp^-1 is applied to the data. The application of using this to form do a V^-1/2 decorrelation is probably medium risk. But this has yet to be proven, and results coming from V^-1/2 should be interpreted with caution.

Note for future: Although the V matrix should be Hermitian by construction, in practice there are precision issues and the Hermiticity is violated at ~ 1 part in 10^15. (Which is ~the expected roundoff error). If something messes up, it may be worth investigating this more.

Note for the future: If this ends up too slow, Cholesky tricks can be employed to speed up the computation by a factor of a few.

Parameters:

• key1, key2 (tuples or lists of tuples) – Tuples containing indices of dataset and baselines for the two input datavectors. If a list of tuples is provided, the baselines in the list will be combined with inverse noise weights.

• exact_norm (boolean, optional) – Exact normalization (see HERA memo #44, Eq. 11 and documentation of q_hat for details).

• pol (str/int/bool, optional) – Polarization parameter to be used for extracting the correct beam. Used only if exact_norm is True.

• model (string, optional) – Type of covariance model to calculate, if not cached. Options=[‘empirical’, ‘dsets’, ‘autos’] How the covariances of the input data should be estimated.

  In ‘dsets’ mode, error bars are estimated from user-provided per baseline and per channel standard deivations.

  If ‘empirical’ is set, then error bars are estimated from the data by averaging the channel-channel covariance of each baseline over time and then applying the appropriate linear transformations to these frequency-domain covariances.

  If ‘autos’ is set, the covariances of the input data over a baseline is estimated from the autocorrelations of the two antennas over channel bandwidth and integration time.

• time_index (int, optional) – Compute covariance at specific time-step. Default: None.

Returns:

V – Bandpower covariance matrix, with dimensions (Ndlys, Ndlys).

Return type:

array_like, complex

iC(key, model='empirical', time_index=None)

Return the inverse covariance matrix, C^-1.

Parameters:

• key (tuple) – Tuple containing indices of dataset and baselines. The first item specifies the index (ID) of a dataset in the collection, while subsequent indices specify the baseline index, in _key2inds format.
• model (string, optional) – Type of covariance model to calculate, if not cached. Options=[‘empirical’, ‘dsets’, ‘autos’] How the covariances of the input data should be estimated. In ‘dsets’ mode, error bars are estimated from user-provided per baseline and per channel standard deivations. If ‘empirical’ is set, then error bars are estimated from the data by averaging the channel-channel covariance of each baseline over time and then applying the appropriate linear transformations to these frequency-domain covariances. If ‘autos’ is set, the covariances of the input data over a baseline is estimated from the autocorrelations of the two antennas over channel bandwidth and integration time.
• time_index (integer, compute covariance at specific time-step)

Returns:

iC – Inverse covariance matrix for specified dataset and baseline.

Return type:

array_like

p_hat(M, q)

Optimal estimate of bandpower p_alpha, defined as p_hat = M q_hat.

Parameters:

• M (array_like) – Normalization matrix, M.
• q (array_like) – Unnormalized bandpowers, hat{q}.

Returns:

p_hat – Optimal estimate of bandpower, hat{p}.

Return type:

array_like

parse_blkey(key)

Parse a dataset + baseline key in the form (dataset, baseline, [pol]) into a dataset index and a baseline(pol) key, where pol is optional.

Parameters:key (tuple) Returns:

• dset (int) – dataset index
• bl (tuple) – baseline (pol) key

pspec(bls1, bls2, dsets, pols, n_dlys=None, input_data_weight='identity', norm='I', taper='none', sampling=False, little_h=True, spw_ranges=None, symmetric_taper=True, baseline_tol=1.0, store_cov=False, store_cov_diag=False, return_q=False, store_window=True, verbose=True, filter_extensions=None, exact_norm=False, history='', r_params=None, cov_model='empirical', known_cov=None)

Estimate the delay power spectrum from a pair of datasets contained in this object, using the optimal quadratic estimator of arXiv:1502.06016.

In this formulation, the power spectrum is proportional to the visibility data via

P = M data_{LH} E data_{RH}

where E contains the data weighting and FT matrices, M is a normalization matrix, and the two separate datasets are denoted as “left-hand” and “right-hand”.

Each power spectrum is generated by taking a baseline (specified by an antenna-pair and polarization key) from bls1 out of dsets[0] and assigning it as data_LH, and a bl from bls2 out of the dsets[1] and assigning it as data_RH.

If the bl chosen from bls1 is (ant1, ant2) and the bl chosen from bls2 is (ant3, ant4), the “baseline-pair” describing their cross multiplication is ((ant1, ant2), (ant3, ant4)).

Parameters:

• bls1, bls2 (list) – List of baseline groups, each group being a list of ant-pair tuples.

• dsets (length-2 tuple or list) – Contains indices of self.dsets to use in forming power spectra, where the first index is for the Left-Hand dataset and second index is used for the Right-Hand dataset (see above).

• pols (tuple or list of tuple) – Contains polarization pairs to use in forming power spectra e.g. (‘XX’,’XX’) or [(‘XX’,’XX’),(‘pI’,’pI’)] or a list of polarization pairs. Individual strings are also supported, and will be expanded into a matching pair of polarizations, e.g. ‘xx’ becomes (‘xx’, ‘xx’).

  If a primary_beam is specified, only equal-polarization pairs can be cross-correlated, as the beam scalar normalization is only implemented in this case. To obtain unnormalized spectra for pairs of different polarizations, set the primary_beam to None.

• n_dlys (list of integer, optional) – The number of delay bins to use. The order in the list corresponds to the order in spw_ranges. Default: None, which then sets n_dlys = number of frequencies.

• input_data_weight (str, optional) – String specifying which weighting matrix to apply to the input data. See the options in the set_R() method for details. Default: ‘identity’.

• norm (str, optional) – String specifying how to choose the normalization matrix, M. See the ‘mode’ argument of get_MW() for options. Default: ‘I’.

• taper (str, optional) – Tapering (window) function to apply to the data. Takes the same arguments as uvtools.dspec.gen_window(). Default: ‘none’.

• sampling (boolean, optional) – Whether output pspec values are samples at various delay bins or are integrated bandpowers over delay bins. Default: False

• little_h (boolean, optional) – Whether to have cosmological length units be h^-1 Mpc or Mpc Default: h^-1 Mpc

• spw_ranges (list of tuples, optional) – A list of spectral window channel ranges to select within the total bandwidth of the datasets, each of which forms an independent power spectrum estimate. Example: [(220, 320), (650, 775)]. Each tuple should contain a start (inclusive) and stop (exclusive) channel used to index the freq_array of each dataset. The default (None) is to use the entire band provided in each dataset.

• symmetric_taper (bool, optional) – Specify if taper should be applied symmetrically to K-matrix (if true) or on the left (if False). Default: True

• baseline_tol (float, optional) – Distance tolerance for notion of baseline “redundancy” in meters. Default: 1.0.

• store_cov (bool, optional) – If True, calculate an analytic covariance between bandpowers given an input visibility noise model, and store the output in the UVPSpec object.

• store_cov_diag (bool, optional) – If True, store the square root of the diagonal of the output covariance matrix calculated by using get_analytic_covariance(). The error bars will be stored in the form of: sqrt(diag(cov_array_real)) + 1.j*sqrt(diag(cov_array_imag)). It’s a way to save the disk space since the whole cov_array data with a size of Ndlys x Ndlys x Ntimes x Nblpairs x Nspws is too large.

• return_q (bool, optional) – If True, return the results (delay spectra and covariance matrices) for the unnormalized bandpowers in the UVPSpec object.

• store_window (bool, optional) – If True, store the window function of the bandpowers. Default: True

• cov_model (string, optional) – Type of covariance model to calculate, if not cached.

  Options=[‘empirical’, ‘dsets’, ‘autos’, ‘foreground_dependent’, (other model names in known_cov)]

  In ‘dsets’ mode, error bars are estimated from user-provided per baseline and per channel standard deivations.

  In ‘empirical’ mode, error bars are estimated from the data by averaging the channel-channel covariance of each baseline over time and then applying the appropriate linear transformations to these frequency-domain covariances.

  In ‘autos’ mode, the covariances of the input data over a baseline is estimated from the autocorrelations of the two antennas forming the baseline across channel bandwidth and integration time.

  In ‘foreground_dependent’ mode, it involves using auto-correlation amplitudes to model the input noise covariance and visibility outer products to model the input systematics covariance.

  For more details see ds.get_analytic_covariance().

• known_cov (dicts of input covariance matrices) – known_cov has a type {Ckey:covariance}, which is the same as ds._C. The matrices stored in known_cov are constructed outside the PSpecData object, unlike those in ds._C which are constructed internally.

  The Ckey should conform to: (dset_pair_index, blpair_int, model, time_index, conj_1, conj_2), e.g. ((0, 1), ((25,37,”xx”), (25, 37, “xx”)), ‘empirical’, False, True), while covariance are ndarrays with shape (Nfreqs, Nfreqs).

  Also see PSpecData.set_C() for more details.

• verbose (bool, optional) – If True, print progress, warnings and debugging info to stdout.

• filter_extensions (list of 2-tuple or 2-list, optional) – Set number of channels to extend filtering width.

• exact_norm (bool, optional) – If True, estimates power spectrum using Q instead of Q_alt (HERA memo #44). The default options is False. Beware that computing a power spectrum when exact_norm is set to False runs two times faster than setting it to True.

• history (str, optional) – History string to attach to UVPSpec object

• r_params (dictionary with parameters for weighting matrix.) – Proper fields and formats depend on the mode of data_weighting.

  For data_weighting set to ‘dayenu’, r_params should be a dict with the following fields:

  filter_centers, a list of floats (or float) specifying the (delay) channel numbers at which to center filtering windows. Can specify fractional channel number.

  filter_half_widths, a list of floats (or float) specifying the width of each filter window in (delay) channel numbers. Can specify fractional channel number.

  filter_factors, list of floats (or float) specifying how much power within each filter window is to be suppressed.

  Absence of an r_params dictionary will result in an error.

Returns:

uvp – Instance of UVPSpec that holds the normalized output power spectrum data.

Return type:

UVPSpec object

Examples

Example 1: No grouping; i.e. each baseline is its own group, no brackets needed for each bl. If:

A = (1, 2); B = (2, 3); C = (3, 4); D = (4, 5); E = (5, 6); F = (6, 7)

and:

bls1 = [ A, B, C ]
bls2 = [ D, E, F ]

then:

blpairs = [ (A, D), (B, E), (C, F) ]

Example 2: Grouping; blpairs come in lists of blgroups, which are considered “grouped” in OQE. If:

bls1 = [ [A, B], [C, D] ]
bls2 = [ [C, D], [E, F] ]

then:

blpairs = [ [(A, C), (B, D)], [(C, E), (D, F)] ]

Example 3: Mixed grouping; i.e. some blpairs are grouped, others are not. If:

bls1 = [ [A, B], C ]
bls2 = [ [D, E], F ]

then:

blpairs = [ [(A, D), (B, E)], (C, F)]

q_hat(key1, key2, allow_fft=False, exact_norm=False, pol=False)

If exact_norm is False: Construct an unnormalized bandpower, q_hat, from a given pair of visibility vectors. Returns the following quantity:

hat{q}_a = (1/2) conj(x_1) R_1 Q^alt_a R_2 x_2

Note that the R matrix need not be set to C^-1. This is something that is set by the user in the set_R method.

This is related to Equation 13 of arXiv:1502.06016. However, notice that there is a Q^alt_a instead of Q_a. The latter is defined as Q_a equiv dC/dp_a. Since this is the derivative of the covariance with respect to the power spectrum, it contains factors of the primary beam etc. Q^alt_a strips away all of this, leaving only the barebones job of taking a Fourier transform. See HERA memo #44 for details.

This function uses the state of self.data_weighting and self.taper in constructing q_hat. See PSpecData.pspec for details.

If exact_norm is True: Takes beam factors into account (Eq. 14 in HERA memo #44)

Parameters:

• key1, key2 (tuples or lists of tuples) – Tuples containing indices of dataset and baselines for the two input datavectors for each power spectrum estimate. q_a formed from key1, key2 If a list of tuples is provided, the baselines in the list will be combined with inverse noise weights.
• allow_fft (bool, optional) – Whether to use a fast FFT summation trick to construct q_hat, or a simpler brute-force matrix multiplication. The FFT method assumes a delta-fn bin in delay space. It also only works if the number of delay bins is equal to the number of frequencies. Default: False.
• exact_norm (bool, optional) – If True, beam and spectral window factors are taken in the computation of Q_matrix (dC/dp = Q, and not Q_alt) (HERA memo #44, Eq. 11). Q matrix, for each delay mode, is weighted by the integral of beam over theta,phi.
• pol (str/int/bool, optional) – Used only if exact_norm is True. This argument is passed to get_integral_beam to extract the requested beam polarization. Default is the first polarization passed to pspec.

Returns:

q_hat – Unnormalized/normalized bandpowers

Return type:

array_like

rephase_to_dset(dset_index=0, inplace=True)

Rephase visibility data in self.dsets to the LST grid of dset[dset_index] using hera_cal.utils.lst_rephase.

Each integration in all other dsets is phased to the center of the corresponding LST bin (by index) in dset[dset_index].

Will only phase if the dataset’s phase type is ‘drift’. This is because the rephasing algorithm assumes the data is drift-phased when applying the phasor term.

Note that PSpecData.Jy_to_mK() must be run after rephase_to_dset(), if one intends to use the former capability at any point.

Parameters:

• dset_index (int or str) – Index or label of dataset in self.dset to phase other datasets to.
• inplace (bool, optional) – If True, edits data in dsets in-memory. Else, makes a copy of dsets, edits data in the copy and returns to user.

Returns:

• if inplace – return new_dsets
• else – return None

scalar(polpair, little_h=True, num_steps=2000, beam=None, taper_override='no_override', exact_norm=False)

Computes the scalar function to convert a power spectrum estimate in “telescope units” to cosmological units, using self.spw_range to set spectral window.

See arxiv:1304.4991 and HERA memo #27 for details.

This function uses the state of self.taper in constructing scalar. See PSpecData.pspec for details.

Parameters:

• polpair (tuple, int, or str) – Which pair of polarizations to compute the beam scalar for, e.g. (‘pI’, ‘pI’) or (‘XX’, ‘YY’). If string, will assume that the specified polarization is to be cross-correlated with itself, e.g. ‘XX’ implies (‘XX’, ‘XX’).
• little_h (boolean, optional) – Whether to have cosmological length units be h^-1 Mpc or Mpc Default: h^-1 Mpc
• num_steps (int, optional) – Number of steps to use when interpolating primary beams for numerical integral Default: 10000
• beam (PSpecBeam object) – Option to use a manually-fed PSpecBeam object instead of using self.primary_beam.
• taper_override (str, optional) – Option to override the taper chosen in self.taper (does not overwrite self.taper; just applies to this function). Default: no_override
• exact_norm (boolean, optional) – If True, scalar would just be the X2Y term, as the beam and spectral terms are taken into account while constructing power spectrum.

Returns:

scalar – [int dnu (Omega_PP / Omega_P^2) ( B_PP / B_P^2 ) / (X^2 Y)]^-1 in h^-3 Mpc^3 or Mpc^3.

Return type:

float

scalar_delay_adjustment(key1=None, key2=None, sampling=False, Gv=None, Hv=None)

Computes an adjustment factor for the pspec scalar that is needed when the number of delay bins is not equal to the number of frequency channels.

This adjustment is necessary because sum_gamma tr[Q^alt_alpha Q^alt_gamma] = N_freq**2 is something that is true only when N_freqs = N_dlys.

If the data weighting is equal to “identity”, then the result is independent of alpha, but is no longer given by N_freq**2. (Nor is it just N_dlys**2!)

If the data weighting is not equal to “identity” then we generally need a separate scalar adjustment for each alpha.

This function uses the state of self.taper in constructing adjustment. See PSpecData.pspec for details.

Parameters:

• key1, key2 (tuples or lists of tuples, optional) – Tuples containing indices of dataset and baselines for the two input datavectors. If a list of tuples is provided, the baselines in the list will be combined with inverse noise weights. If Gv and Hv are specified, these arguments will be ignored. Default: None.
• sampling (boolean, optional) – Whether to sample the power spectrum or to assume integrated bands over wide delay bins. Default: False
• Gv, Hv (array_like, optional) – If specified, use these arrays instead of calling self.get_G() and self.get_H(). Using precomputed Gv and Hv will speed up this function significantly. Default: None.

Returns:

adjustment – 1d array of floats with length spw_Ndlys otherwise.

Return type:

float if the data_weighting is ‘identity’

set_C(cov)

Set the cached covariance matrix to a set of user-provided values.

Parameters:cov (dict) – Covariance keys and ndarrays. The key should conform to (dset_pair_index, blpair_int, model, time_index, conj_1, conj_2). e.g. ((0, 1), ((25,37,”xx”), (25, 37, “xx”)), ‘empirical’, False, True) while the ndarrays should have shape (spw_Nfreqs, spw_Nfreqs)

set_Ndlys(ndlys=None)

Set the number of delay bins used.

Parameters:ndlys (integer) – Number of delay bins. Default: None, sets number of delay bins equal to the number of frequency channels in the current spw

set_R(d)

Set the data-weighting matrix for a given dataset and baseline to a specified value for later use in q_hat.

Parameters:

d (dict) – Dictionary containing data to insert into data-weighting R matrix cache. Keys are tuples with the following form key = (dset_index, bl_ant_pair_pol_tuple, data_weighting, taper) Example: (0, (37, 38, ‘xx’), ‘bh’)

If data_weight == ‘dayenu’ then additional elements are appended: key + (filter_extension, spw_Nfreqs, symmetric_taper)

set_filter_extension(filter_extension)

Set extensions to filtering matrix

Parameters:filter_extension (2-tuple or 2-list) – must be integers. Specify how many channels below spw_min/max filter will be applied to data. filter_extensions will be clipped to not extend beyond data range.

set_iC(d)

Set the cached inverse covariance matrix for a given dataset and baseline to a specified value. For now, you should already have applied weights to this matrix.

Parameters:d (dict) – Dictionary containing data to insert into inverse covariance matrix cache. Keys are tuples, following the same format as the input to self.iC().

set_r_param(key, r_params)

Set the weighting parameters for baseline at (dset,bl, [pol])

Parameters:

• key (tuple) – Key in the format: (dset, bl, [pol]) dset is the index of the dataset, bl is a 2-tuple, pol is a float or string specifying polarization.

• r_params (dict) – Dict containing parameters for weighting matrix. Proper fields and formats depend on the mode of data_weighting.

  For data_weighting set to dayenu, this is a dictionary with the following fields:

  filter_centers, list of floats (or float) specifying the (delay) channel numbers at which to center filtering windows. Can specify fractional channel number.

  filter_half_widths, list of floats (or float) specifying the width of each filter window in (delay) channel numbers. Can specify fractional channel number.

  filter_factors, list of floats (or float) specifying how much power within each filter window is to be suppressed.

  Absence of an r_params dictionary will result in an error.

set_spw(spw_range, ndlys=None)

Set the spectral window range.

Parameters:

• spw_range (tuple, contains start and end of spw in channel indices) – used to slice the frequency array
• ndlys (integer) – Number of delay bins. Default: None, sets number of delay bins equal to the number of frequency channels in the spw.

set_symmetric_taper(use_symmetric_taper)

Set the symmetric taper parameter If true, square R matrix will be computed as R=sqrtT K sqrt T where sqrtT is a diagonal matrix with the square root of the taper. This is only possible when K is a square-matrix (no filter extensions).

If set to false, then the R-matrix will implement the taper as R = sqrtT ** 2 K :Parameters: use_taper (bool,) – do you want to use a symmetric taper? True or False?

set_taper(taper)

Set data tapering type.

Parameters:taper (str) – Type of data tapering. See uvtools.dspec.gen_window for options.

set_weighting(data_weighting)

Set data weighting type.

Parameters:data_weighting (str) – Type of data weightings. Options=[‘identity’, ‘iC’, dayenu]

trim_dset_lsts(lst_tol=6)

Assuming all datasets in self.dsets are locked to the same LST grid (but each may have a constant offset), trim LSTs from each dset that aren’t found in all other dsets (within some decimal tolerance specified by lst_tol).

Warning: this edits the data in dsets in-place, and is not reversible.

Parameters:lst_tol (float) – Decimal tolerance [radians] for comparing float-valued LST bins.

units(little_h=True)

Return the units of the power spectrum. These are inferred from the units reported by the input visibilities (UVData objects).

Parameters:little_h (boolean, optional) – Whether to have cosmological length units be h^-1 Mpc or Mpc Default: h^-1 Mpc Returns:pspec_units – Units of the power spectrum that is returned by pspec(). Return type:str

validate_datasets(verbose=True)

Validate stored datasets and weights to make sure they are consistent with one another (e.g. have the same shape, baselines etc.).

validate_pol(dsets, pol_pair)

Validate polarization and returns the index of the datasets so that the polarization pair is consistent with the UVData objects.

Parameters:

• dsets (length-2 list or length-2 tuple of integers or str) – Contains indices of self.dsets to use in forming power spectra, where the first index is for the Left-Hand dataset and second index is used for the Right-Hand dataset (see above).
• pol_pair (length-2 tuple of integers or strings) – Contains polarization pair which will be used in estiamting the power spectrum e,g (-5, -5) or (‘xy’, ‘xy’). Only auto-polarization pairs are implemented for the time being.

Returns:

valid – True if the UVData objects polarizations are consistent with the pol_pair (user specified polarizations) else False.

Return type:

boolean

w(key, include_extension=False)

Get weights for a given dataset and baseline, as specified in a standard key format.

Parameters:

• key (tuple) – Tuple containing dataset ID and baseline index. The first element of the tuple is the dataset index, and the subsequent elements are the baseline ID.
• include_extension (bool (optional)) – default=False If True, extend spw to include filtering window extensions.

Returns:

w – Array of weights for the requested UVData dataset and baseline.

Return type:

array_like

x(key, include_extension=False)

Get data for a given dataset and baseline, as specified in a standard key format.

Parameters:

• key (tuple) – Tuple containing dataset ID and baseline index. The first element of the tuple is the dataset index (or label), and the subsequent elements are the baseline ID.
• include_extension (bool (optional)) – default=False If True, extend spw to include filtering window extensions.

Returns:

x – Array of data from the requested UVData dataset and baseline.

Return type:

array_like

UVPSpec Class

UVPSpec objects contain power spectra and associated metadata.

class hera_pspec.UVPSpec

An object for storing power spectra generated by hera_pspec and a file-format for its data and meta-data.

antnums_to_bl(antnums)

Convert tuple of antenna numbers to baseline integer.

Parameters:antnums (tuple) – Tuple containing integer antenna numbers for a baseline. Ex. (ant1, ant2) Returns:bl – Baseline integer. Return type:i6 integer

antnums_to_blpair(antnums)

Convert nested tuple of antenna numbers to baseline-pair integer.

Parameters:antnums (tuple) – nested tuple containing integer antenna numbers for a baseline-pair. Ex. ((ant1, ant2), (ant3, ant4)) Returns:blpair – baseline-pair integer Return type:i12 integer

average_spectra(blpair_groups=None, time_avg=False, blpair_weights=None, error_field=None, error_weights=None, inplace=True, add_to_history='')

Average power spectra across the baseline-pair-time axis, weighted by each spectrum’s integration time.

This is an “incoherent” average, in the sense that this averages power spectra, rather than visibility data. The nsample_array and integration_array will be updated to reflect the averaging.

In the case of averaging across baseline pairs, the resultant averaged spectrum is assigned to the zeroth blpair in the group. In the case of time averaging, the time and LST arrays are assigned to the mean of the averaging window.

Note that this is designed to be separate from spherical binning in k: here we are not connecting k_perp modes to k_para modes. However, if blpairs holds groups of iso baseline separation, then this is equivalent to cylindrical binning in 3D k-space.

If you want help constructing baseline-pair groups from baseline groups, see self.get_blpair_groups_from_bl_groups.

Parameters:

blpair_groups (list) – List of list of baseline-pair group tuples or integers. All power spectra in a baseline-pair group are averaged together. If a baseline-pair exists in more than one group, a warning is raised.

Examples:

blpair_groups = [ [((1, 2), (1, 2)), ((2, 3), (2, 3))],
                  [((4, 6), (4, 6))]]

blpair_groups = [ [1002001002, 2003002003], [4006004006] ]

time_avg : bool, optional

If True, average power spectra across the time axis. Default: False.

blpair_weights : list, optional

List of float or int weights dictating the relative weight of each baseline-pair when performing the average.

This is useful for bootstrapping. This should have the same shape as blpair_groups if specified. The weights are automatically normalized within each baseline-pair group. Default: None (all baseline pairs have unity weights).

error_field: string or list, optional

If errorbars have been entered into stats_array, will do a weighted sum to shrink the error bars down to the size of the averaged data_array. Error_field strings be keys of stats_array. If list, does this for every specified key. Every stats_array key that is not specified is thrown out of the new averaged object.

error_weights: string, optional

error_weights specify which kind of errors we use for weights during averaging power spectra. The weights are defined as $w_i = 1/ sigma_i^2$, where $sigma_i$ is taken from the relevant field of stats_array.

If error_weight is set to None, which means we just use the integration time as weights. If error_weights is specified, then it also gets appended to error_field as a list. Default: None.

inplace : bool, optional

If True, edit data in self, else make a copy and return. Default: True.

add_to_history : str, optional

Added text to add to file history.

Notes

Currently, every baseline-pair in a blpair group must have the same Ntimes, unless time_avg=True. Future versions may support baseline-pair averaging of heterogeneous time arrays. This includes the scenario of repeated blpairs (e.g. in bootstrapping), which will return multiple copies of their time_array.

bl_to_antnums(bl)

Convert baseline (antenna-pair) integer to nested tuple of antenna numbers.

Parameters:bl (i6 int) – baseline integer ID Returns:antnums – tuple containing baseline antenna numbers. Ex. (ant1, ant2) Return type:tuple

blpair_to_antnums(blpair)

Convert baseline-pair integer to nested tuple of antenna numbers.

Parameters:blpair (i12 int) – baseline-pair integer ID Returns:antnums – Nested tuple containing baseline-pair antenna numbers. Ex. ((ant1, ant2), (ant3, ant4)) Return type:tuple

blpair_to_indices(blpair)

Convert a baseline-pair nested tuple ((ant1, ant2), (ant3, ant4)) or a baseline-pair integer into indices to index the blpairts axis of data_array.

Parameters:blpair (tuples or ints) – Nested tuple or blpair i12 integer, Ex. ((1, 2), (3, 4)), or list of blpairs.

check(just_meta=False)

Run checks to make sure metadata and data arrays are properly defined and in the right format. Raises AssertionErrors if checks fail.

Parameters:just_meta (bool, optional) – If True, only check metadata. Default: False.

compute_scalar(spw, polpair, num_steps=1000, little_h=True, noise_scalar=False)

Compute power spectrum normalization scalar given an adopted cosmology and a beam model. See pspecbeam.PSpecBeamBase.compute_pspec_scalar for details.

Parameters:

• spw (integer) – Spectral window selection.
• polpair (tuple or int or str) – Which polarization pair to calculate the scalar for.
• num_steps (int, optional) – Number of integration bins along frequency in computing scalar. Default: 1000.
• little_h (bool, optional) – Whether to use h^-1 units or not. Default: True.
• noise_scalar (bool, optional) – If True calculate noise pspec scalar, else calculate normal pspec scalar. See pspecbeam.py for difference between normal scalar and noise scalar. Default: False.

Returns:

scalar – Power spectrum normalization scalar.

Return type:

float

convert_to_deltasq(inplace=True)

Convert from P(k) to Delta^2(k) by multiplying by k^3 / (2pi^2).

The units of the output is therefore the current units (self.units) times (h^3) Mpc^-3, where the h^3 is only included if h^-3 is in self.norm_units

Parameters:

• inplace (bool, optional) – If True edit and overwrite arrays in self, else make a copy of self and return. Default: True.
• Notes
• ——
• Alters data_array, stats_array, and cov_array if present,
• as well as metadata like norm_units

fold_spectra()

Average bandpowers from matching positive and negative delay bins onto a purely positive delay axis. Negative delay bins are still populated, but are filled with zero power. This is an in-place operation.

Will only work if self.folded == False, i.e. data is currently unfolded across negative and positive delay. Because this averages the data, the nsample array is multiplied by a factor of 2.

WARNING: This operation cannot be undone.

generate_noise_spectra(spw, polpair, Tsys, blpairs=None, little_h=True, form='Pk', num_steps=2000, component='real')

Generate the expected RMS noise power spectrum given a selection of spectral window, system temp. [K], and polarization. This estimate is constructed as:

P_N = scalar * (Tsys * 1e3)^2 / (integration_time) / sqrt(Nincoherent)

[mK^2 h^-3 Mpc^3]

where scalar is the cosmological and beam scalar (i.e. X2Y * Omega_eff) calculated from pspecbeam with noise_scalar = True, integration_time is in seconds and comes from self.integration_array and Nincoherent is the number of incoherent averaging samples and comes from self.nsample_array. If component is real or imag, P_N is divided by an additional factor of sqrt(2).

If the polarizations specified are pseudo Stokes pol (I, Q, U or V) then an extra factor of 2 is divided.

If form is DelSq then a factor of |k|^3 / (2pi^2) is multiplied.

If real is True, a factor of sqrt(2) is divided to account for discarding imaginary noise component.

For more details, see hera_pspec.noise.Sensitivity.calc_P_N, and Cheng et al. 2018.

The default units of P_N are mK^2 h^-3 Mpc^3 (if little_h=True and form=’Pk’), but is different if those kwargs are altered.

Parameters:

• spw (int) – Spectral window index to generate noise curve for.
• polpair (tuple or int or str) – Polarization-pair selection in form of tuple (e.g. (‘I’, ‘Q’)) or polpair int. Strings are expanded into polarization pairs, e.g. ‘XX’ becomes (‘XX,’XX’).
• Tsys (dictionary, float or array) – System temperature in Kelvin for each blpair. Key is blpair-integer, value is Tsys float or ndarray. If fed as an ndarray, shape=(Ntimes,)
• blpairs (list) – List of unique blair tuples or i12 integers to calculate noise spectrum for default is to calculate for baseline pairs.
• little_h (bool, optional) – Whether to have cosmological length units be h^-1 Mpc or Mpc Default: h^-1 Mpc
• form (str, optional) – Form of power spectra, whetehr P(k) or Delta^2(k). Valid options are ‘Pk’ or ‘DelSq’. Default: ‘Pk’.
• num_steps (int, optional) – Number of frequency bins to use in integrating power spectrum scalar in pspecbeam. Default: 2000.
• component (str) – Options=[‘real’, ‘imag’, ‘abs’]. If component is real or imag, divide by an extra factor of sqrt(2).

Returns:

P_N – Dictionary containing blpair integers as keys and float ndarrays of noise power spectra as values, with ndarrays having shape (Ntimes, Ndlys).

Return type:

dict

get_ENU_bl_vecs()

Return baseline vector array in TOPO (ENU) frame in meters, with matched ordering of self.bl_vecs.

Returns:blvecs – Array of baseline vectors. Return type:array_like

get_all_keys()

Returns a list of all possible tuple keys in the data_array, in the format:

(spectral window, baseline-pair, polarization-string)

get_blpair_blvecs(use_second_bl=False)

For each baseline-pair, get the baseline vector in ENU frame of the first baseline in the pair. Ordering matches self.get_blpairs()

Parameters:use_second_bl (bool) – If True, return the vector of the second bl in the blpair. Returns:vector in ENU frame of the baseline Return type:ndarray shape=(Nblpairs,)

get_blpair_groups_from_bl_groups(blgroups, only_pairs_in_bls=False)

Get baseline pair matches from a list of baseline groups.

Parameters:

• blgroups (list) – List of baseline groups, which themselves are lists of baseline tuples or baseline i6 integers. Ex: [ [(1, 2), (2, 3), (3, 4)], [(1, 4), (2, 5)] ]
• only_pairs_in_bls (bool, optional) – If True, select only baseline-pairs whose first _and_ second baseline are both found in each baseline group. Default: False.

Returns:

blpair_groups – List of blpair groups, which themselves are lists of blpair integers.

Return type:

list

get_blpair_seps()

For each baseline-pair, get the average baseline separation in ENU frame in meters between its bl1 and bl2. Ordering matches self.get_blpairs()

Returns:blp_avg_sep – shape=(Nblpairts,) Return type:float ndarray

get_blpairs()

Returns a list of all blpair tuples in the data_array.

get_cov(key, component='real', omit_flags=False)

Slice into covariance array with a specified data key in the format (spw, ((ant1, ant2),(ant3, ant4)), (pol1, pol2))

or

(spw, blpair-integer, pol-integer)

where spw is the spectral window integer, ant1 etc. are integers, and pol is either a tuple of polarization strings (ex. (‘XX’, ‘YY’)) or integers (ex. (-5,-5)).

Parameters:

• key (tuple) – Contains the baseline-pair, spw, polpair keys
• component (str) – “real” or “imag”. Indicating which cov_array the function calls.
• omit_flags (bool, optional) – If True, remove time integrations (or spectra) that came from visibility data that were completely flagged across the spectral window (i.e. integration == 0).

Returns:

data – Shape (Ntimes, Ndlys, Ndlys)

Return type:

complex ndarray

get_data(key, omit_flags=False)

Slice into data_array with a specified data key in the format

(spw, ((ant1, ant2), (ant3, ant4)), (pol1, pol2))

or

(spw, blpair-integer, polpair)

where spw is the spectral window integer, ant1 etc. are integers, and polpair is either a tuple of polarization strings/ints or a polarizatio-pair integer.

Parameters:

• key (tuple) – Contains the baseline-pair, spw, polpair keys
• omit_flags (bool, optional) – If True, remove time integrations (or spectra) that came from visibility data that were completely flagged across the spectral window (i.e. integration == 0).

Returns:

data – Shape (Ntimes, Ndlys)

Return type:

complex ndarray

get_dlys(key)

Get array of delays given a spectral window selection.

Parameters:key (int, or tuple with integer) – Spectral window selection Returns:dlys Return type:float ndarray, contains delays of pspectra in nanosec, given spw

get_integrations(key, omit_flags=False)

Slice into integration_array with a specified data key in the format:

(spw, ((ant1, ant2), (ant3, ant4)), (pol1, pol2))

or

(spw, blpair-integer, polpair-integer)

where spw is the spectral window integer, ant1 etc. are integers, and polpair is either a polarization pair tuple or integer.

Parameters:

• key (tuple) – Contains the baseline-pair, spw, polpair keys
• omit_flags (bool, optional) – If True, remove time integrations (or spectra) that came from visibility data that were completely flagged across the spectral window (i.e. integration == 0).

Returns:

data – Has shape (Ntimes,)

Return type:

float ndarray

get_kparas(spw, little_h=True)

Get radial (parallel) cosmological wavevectors for power spectra given an adopted cosmology and a spw selection.

Parameters:

• spw (int) – Choice of spectral window
• little_h (boolean, optional) – Whether to have cosmological length units be h^-1 Mpc or Mpc Default: h^-1 Mpc
• Returns (k_perp, k_para)
• ——-
• k_para (float ndarray) – Radial wave-vectors, shape=(Ndlys given spw)

get_kperps(spw, little_h=True)

Get transverse (perpendicular) cosmological wavevector for each baseline-pair given an adopted cosmology and a spw selection.

Parameters:

• spw (int, choice of spectral window)
• little_h (boolean, optional) – Whether to have cosmological length units be h^-1 Mpc or Mpc Default: h^-1 Mpc

Returns:

k_perp

Return type:

float ndarray, transverse wave-vectors, shape=(Nblpairs,)

get_nsamples(key, omit_flags=False)

Slice into nsample_array with a specified data key in the format

(spw, ((ant1, ant2), (ant3, ant4)), (pol1, pol2))

or

(spw, blpair-integer, polpair-integer)

where spw is the spectral window integer, ant1 etc. are integers, and polpair is either a polarization pair tuple or integer.

Parameters:

• key (tuple) – Contains the baseline-pair, spw, polpair keys
• omit_flags (bool, optional) – If True, remove time integrations (or spectra) that came from visibility data that were completely flagged across the spectral window (i.e. integration == 0).

Returns:

data

Return type:

float ndarray with shape (Ntimes,)

get_polpairs()

Returns a list of all unique polarization pairs in the data_array.

get_r_params()

Return an r_params dictionary.

In a PSpecData object, the r_params are stored as a dictionary with one entry per baseline.

In a UVPSpec object, the dictionary is compressed so that a single r_param entry correspondsto multiple baselines and is stored as a JSON format string.

This function reads the compressed string and returns the dictionary with the correct following fields and structure.

Returns:r_params – Dictionary of r_params for this object. Return type:dict

get_spw_ranges(spws=None)

Get the frequency range and Nfreqs of each spectral window in the object.

Parameters:spws (list of spw integers, optional) – Default is all spws present in data. Returns:spw_ranges – Tuples: (freq_start, freq_end, Nfreqs, Ndlys) in Hz. Contains start, stop and Nfreqs and Ndlys [Hz] of each spectral window. To turn this into the frequency array of the spectral window use spw_freqs = np.linspace(*spw_range[:3], endpoint=False) Return type:list of len-4 tuples

get_stats(stat, key, omit_flags=False)

Returns a statistic from the stats_array dictionary.

Parameters:

• stat (str) – The statistic to return.
• key (tuple) – A spw-blpair-polpair key parseable by UVPSpec.key_to_indices.
• omit_flags (bool, optional) – If True, remove time integrations (or spectra) that came from visibility data that were completely flagged across the spectral window (i.e. integration == 0).

Returns:

statistic – A 2D (Ntimes, Ndlys) complex array holding desired bandpower statistic.

Return type:

array_like

get_wgts(key, omit_flags=False)

Slice into wgt_array with a specified data key in the format

(spw, ((ant1, ant2), (ant3, ant4)), (pol1, pol2))

or

(spw, blpair-integer, polpair-integer)

where spw is the spectral window integer, ant1 etc. are integers, and polpair is either a polarization pair tuple or integer.

Parameters:

• key (tuple) – Contains the baseline-pair, spw, polpair keys
• omit_flags (bool, optional) – If True, remove time integrations (or spectra) that came from visibility data that were completely flagged across the spectral window (i.e. integration == 0).

Returns:

wgts – Has shape (Ntimes, Nfreqs, 2), where the last axis holds [wgt_1, wgt_2] in that order.

Return type:

float ndarray

get_window_function(key, omit_flags=False)

Slice into window_function array with a specified data key in the format (spw, ((ant1, ant2),(ant3, ant4)), (pol1, pol2))

or

(spw, blpair-integer, pol-integer)

where spw is the spectral window integer, ant1 etc. are integers, and pol is either a tuple of polarization strings (ex. (‘XX’, ‘YY’)) or integers (ex. (-5,-5)).

Parameters:

• key (tuple) – Contains the baseline-pair, spw, polpair keys
• omit_flags (bool, optional) – If True, remove time integrations (or spectra) that came from visibility data that were completely flagged across the spectral window (i.e. integration == 0).

Returns:

data – Shape (Ntimes, Ndlys, Ndlys)

Return type:

complex ndarray

key_to_indices(key, omit_flags=False)

Convert a data key into relevant slice arrays. A data key takes the form

(spw_integer, ((ant1, ant2), (ant3, ant4)), (pol1, pol2))

or

(spw, blpair-integer, pol-integer)

where spw is the spectral window integer, ant1 etc. are integers, and polpair is either a polarization-pair integer or tuple.

The key can also be a dictionary in the form:

key = {
  'spw' : spw_integer,
  'blpair' : ((ant1, ant2), (ant3, ant4))
  'polpair' : (pol1, pol2)
  }

and it will parse the dictionary for you.

Parameters:

• key (tuple) – Baseline-pair, spw, and pol-pair key.
• omit_flags (bool, optional) – If True, remove time integrations (or spectra) that came from visibility data that were completely flagged across the spectral window (i.e. integration == 0).

Returns:

• spw (int) – Spectral window index.
• blpairts (list) – List of integers to apply along blpairts axis.
• polpair (int) – Polarization pair index.

polpair_to_indices(polpair)

Map a polarization-pair integer or tuple to its index in polpair_array.

Parameters:

polpair ((list of) int or tuple or str) – Polarization-pair integer or tuple of polarization strings/ints.

Alternatively, if a single string is given, will assume that the specified polarization is to be cross-correlated withitself, e.g. ‘XX’ implies (‘XX’, ‘XX’).

Returns:

indices – Index of polpair in polpair_array.

Return type:

int

read_from_group(grp, just_meta=False, spws=None, bls=None, blpairs=None, times=None, lsts=None, polpairs=None, only_pairs_in_bls=False)

Clear current UVPSpec object and load in data from specified HDF5 group.

Parameters:

• grp (HDF5 group) – HDF5 group to load data from.
• just_meta (bool, optional) – If True, read-in metadata but ignore data, wgts and integration arrays. Default: False.
• spws (list of tuple, optional) – List of spectral window integers to select. Default: None (loads all channels).
• bls (list of i6 baseline integers or baseline tuples) – Select all baseline-pairs whose first _or_ second baseline are in the list. This changes if only_pairs_in_bls == True. Example tuple: (2, 3). Default: None (loads all bl pairs).
• blpairs (list of baseline-pair tuples or integers) – List of baseline pairs to keep. If bls is also fed, this list is concatenated onto the baseline-pair list constructed from from the bls selection.
• times (float ndarray) – Times from the time_avg_array to keep.
• lsts (float ndarray) – lsts from the lst_avg_array to keep.
• polpairs (list of tuple or int or str) – List of polarization-pair integers, tuples, or strings to keep. Strings are expanded into polarization pairs, and are not treated as individual polarizations, e.g. ‘XX’ becomes (‘XX,’XX’).
• only_pairs_in_bls (bool, optional) – If True, keep only baseline-pairs whose first _and_ second baseline are found in the ‘bls’ list. Default: False.

read_hdf5(filepath, just_meta=False, spws=None, bls=None, blpairs=None, times=None, lsts=None, polpairs=None, only_pairs_in_bls=False)

Clear current UVPSpec object and load in data from an HDF5 file.

Parameters:

• filepath (str) – Path to HDF5 file to read from.

• just_meta (bool, optional) – If True, read-in metadata but ignore data, wgts and integration arrays. Default: False.

• spws (list of tuple, optional) – List of spectral window integers to select. Default: None (loads all channels).

• bls (list of i6 baseline integers or baseline tuples) – Select all baseline-pairs whose first _or_ second baseline are in the list. This changes if only_pairs_in_bls == True. Example tuple: (2, 3). Default: None (loads all bl pairs).

• blpairs (list of baseline-pair tuples or integers) – List of baseline pairs to keep. If bls is also fed, this list is concatenated onto the baseline-pair list constructed from from the bls selection.

• times (float ndarray) – Times from the time_avg_array to keep.

• lsts (float ndarray) – lsts from the lst_avg_array to keep.

• polpairs (list of polpair ints or tuples or str) – List of polarization-pair integers or tuples to keep.

  Strings are expanded into polarization pairs, and are not treated as individual polarizations, e.g. ‘XX’ becomes (‘XX,’XX’).

• only_pairs_in_bls (bool, optional) – If True, keep only baseline-pairs whose first _and_ second baseline are found in the ‘bls’ list. Default: False.

select(spws=None, bls=None, only_pairs_in_bls=True, blpairs=None, times=None, lsts=None, polpairs=None, inplace=True, run_check=True)

Select function for selecting out certain slices of the data.

Parameters:

• spws (list, optional) – List of spectral window integers to select. If None, all will be selected. Default: None.

• bls (list of i6 baseline integers or baseline tuples, optional) – Example: (2,3). Select all baseline-pairs whose first _or_ second baseline are in bls list. This changes if only_pairs_in_bls == True. If None, all baselines are selected (subject to other options). Default: None.

• only_pairs_in_bls (bool, optional) – If True, keep only baseline-pairs whose first _and_ second baseline are found in bls list. Default: True.

• blpairs (list of baseline-pair tuples or integers, optional) – List of baseline-pairs to keep. If bls is also fed, this list is concatenated onto the baseline-pair list constructed from the bls selection. If None, all valid baseline pairs are selected. Default: None.

• times (array_like, optional) – Float ndarray of times from the time_avg_array to select. If None, all times are kept. Default: None.

• lsts (array_like, optional) – Float ndarray of lsts from the lst_avg_array to select. If None, all lsts are kept. Default: None.

• polpairs (list of tuple or int or str, optional) – List of polarization-pairs to keep. If None, all polarizations are kept. Default: None.

  (Strings are expanded into polarization pairs, and are not treated as individual polarizations, e.g. ‘XX’ becomes (‘XX,’XX’).)

• inplace (bool, optional) – If True, edit and overwrite arrays in self, else make a copy of self and return. Default: True.

• run_check (bool, optional) – If True, run uvp.check() after selection

Returns:

uvp – If inplace=False, return a new UVPSpec object containing only the selected data.

Return type:

UVPSpec, optional

set_cosmology(new_cosmo, overwrite=False, new_beam=None, verbose=True)

Set the cosmology for this UVPSpec object by passing an instance of conversions.Cosmo_Conversions and assigning it to self.cosmo, in addition to re-computing power spectrum normalization quantities like self.scalar_array. Because this will attempt to recompute the scalar_array, the beam-related metadata (OmegaP, OmegaPP and beam_freqs) must exist in self.

If they do not, or if you’d like to overwrite them with a new beam, you can pass a UVBeam object or path to a beam file via the new_beam kwarg.

If self.cosmo already exists then you are attempting to overwrite the currently adopted cosmology, which will only proceed if overwrite is True.

If overwrite == True, then this module will overwrite self.cosmo. It will also recompute the power spectrum normalization scalar (using pspecbeam) and will overwrite the values in self.scalar_array. It will then propagate these re-normalization changes into the data_array by multiplying the data by (new_scalar_array / old_scalar_array).

Parameters:

• new_cosmo (Cosmo_Conversions or dict) – Cosmo_Conversions instance or cosmological parameter dictionary. The new cosmology you want to adopt for this UVPSpec object.
• overwrite (bool, optional) – If True, overwrite self.cosmo if it already exists. Default: False.
• new_beam (PSpecBeamBase sublcass or str) – pspecbeam.PSpecBeamUV object or path to beam file. The new beam you want to adopt for this UVPSpec object.
• verbose (bool, optional) – If True, report feedback to stdout. Default: True.

set_stats(stat, key, statistic)

Set a statistic waterfall (Ntimes, Ndlys) in the stats_array given the selection of spw, blpairts and polpair in key.

Parameters:

• stat (string) – Name of the statistic.
• key (tuple) – A spw-blpair-polpair key parseable by UVPSpec.key_to_indices.
• statistic (ndarray) – 2D ndarray with statistics to set. Must conform to shape of the slice of data_array produced by the specified key.

set_stats_slice(stat, m, b, above=False, val=1e+20)

For each baseline, set all delay bins in stats_array that fall above or below y = bl_len * m + b [nanosec] equal to val. Useful for downweighting foregrounds in spherical average.

Parameters:

• stat (str, name of stat in stat_array to set)
• m (float, coefficient of bl_len [nanosec / meter])
• b (float, offset [nanosec])
• above (bool, if True, set stats above line, else set below line)
• val (float, value to replace in stats_array)

spw_indices(spw)

Convert a spectral window integer into a list of indices to index into the spw_array.

Parameters:spw (int or tuple) – Spectral window index, or spw tuple from get_spw_ranges(), or list of either. Returns:spw_indices – Index array for specified spw(s) in spw_array. Return type:array_like

spw_to_dly_indices(spw)

Convert a spectral window integer into a list of indices to index into the spw_dly_array or dly_array.

Parameters:spw (int or tuple) – Spectral window index, or spw tuple from get_spw_ranges(), or list of either. Returns:dly_indices – Index array for specified spw(s) in spw_dly_array. Return type:array_like

spw_to_freq_indices(spw)

Convert a spectral window integer into a list of indices to index into the spw_freq_array or freq_array.

Parameters:spw (int or tuple) – Spectral window index, or spw tuple from get_spw_ranges(), or list of either. Returns:freq_indices – Index array for specified spw(s) in spw_freq_array. Return type:array_like

time_to_indices(time, blpairs=None)

Convert a time [Julian Date] from self.time_avg_array to an array that indexes the elements at which it appears in that array. Can optionally take a blpair selection to further select the indices at which both that time and blpair are satisfied.

Parameters:

• time (float) – Julian Date time from self.time_avg_array, Ex: 2458042.12242
• blpairs (tuple or int, optional) – List of blpair tuples or integers that further selects the elements at which both time and blpairs are satisfied.

Returns:

indices – Contains indices at which selection is satisfied along the blpairts axis.

Return type:

integer ndarray

units

Return power spectrum units. See self.vis_units and self.norm_units.

write_hdf5(filepath, overwrite=False, run_check=True)

Write a UVPSpec object to HDF5 file.

Parameters:

• filepath (str) – Filepath for output file.
• overwrite (bool, optional) – Whether to overwrite output file if it exists. Default: False.
• run_check (bool, optional) – Run UVPSpec validity check before writing to file. Default: True.

write_to_group(group, run_check=True)

Write UVPSpec data into an HDF5 group.

Parameters:

• group (HDF5 group) – The handle of the HDF5 group that the UVPSpec object should be written into.
• run_check (bool, optional) – Whether to run a validity check on the UVPSpec object before writing it to the HDF5 group. Default: True.

PSpecContainer Class

PSpecContainer is a container for organizing collections of UVPSpec objects. It is based on HDF5.

class hera_pspec.PSpecContainer(filename, mode='r', keep_open=True, swmr=False, tsleep=0.1, maxiter=1)

Container class for managing multiple UVPSpec objects.

get_pspec(group, psname=None, **kwargs)

Get a UVPSpec power spectrum object from a given group.

Parameters:

• group (str) – Which group the power spectrum belongs to.
• psname (str, optional) – The name of the power spectrum to return. If None, extract all available power spectra.
• kwargs (dict) – UVPSpec.read_from_group partial IO keyword arguments

Returns:

uvp – The specified power spectrum as a UVPSpec object (or a list of all power spectra in the group, if psname was not specified).

Return type:

UVPSpec or list of UVPSpec

groups()

Return list of groups in the container.

Returns:group_list – List of group names. Return type:list of str

save()

Force HDF5 file to flush to disk.

set_pspec(group, psname, pspec, overwrite=False)

Store a delay power spectrum in the container.

Parameters:

• group (str) – Which group the power spectrum belongs to.
• psname (str or list of str) – The name(s) of the power spectrum to return from within the group.
• pspec (UVPSpec or list of UVPSpec) – Power spectrum object(s) to store in the container.
• overwrite (bool, optional) – If the power spectrum already exists in the file, whether it should overwrite it or raise an error. Default: False (does not overwrite).

spectra(group)

Return list of available power spectra.

Parameters:group (str) – Which group to list power spectra from. Returns:ps_list – List of names of power spectra in the group. Return type:list of str

tree(return_str=True)

Output a string containing a tree diagram of groups and the power spectra that they contain.

Parameters:return_str (bool, optional) – If True, return the tree as a string, otherwise return as a dictionary Returns:Tree structure of HDF5 file Return type:str or dict

Simple plotting functions

The hera_pspec.plot module contains functions for making simple plots of delay power spectra.

The following example plots the power spectra from a UVPSpec object, averaged over baseline-pairs and times.

# Load or generate a UVPSpec object containing delay power spectra
uvp = ...

# Set which baseline-pairs should be included in the plot
blpairs = list(uvp.blpair_array) # This includes all blpairs!

# Plot the delay spectrum, averaged over all blpairs and times
# (for the spectral window with index=0, and polarization 'xx')
ax = hp.plot.delay_spectrum(uvp, [blpairs,], spw=0, pol='xx',
                            average_blpairs=True, average_times=True,
                            delay=False)

# Setting delay=False plots the power spectrum in cosmological units

For a more extensive worked example, see this example Jupyter notebook.

The only plotting function currently available in the hera_pspec.plot module is delay_spectrum().

hera_pspec.plot.delay_spectrum(uvp, blpairs, spw, pol, average_blpairs=False, average_times=False, fold=False, plot_noise=False, delay=True, deltasq=False, legend=False, ax=None, component='real', lines=True, markers=False, error=None, times=None, logscale=True, force_plot=False, label_type='blpairt', plot_stats=None, **kwargs)

Plot a 1D delay spectrum (or spectra) for a group of baselines.

Parameters:

• uvp (UVPspec) – UVPSpec object, containing delay spectra for a set of baseline-pairs, times, polarizations, and spectral windows.

• blpairs (list of tuples or lists of tuples) – List of baseline-pair tuples, or groups of baseline-pair tuples.

• spw, pol (int or str) – Which spectral window and polarization to plot.

• average_blpairs (bool, optional) – If True, average over the baseline pairs within each group.

• average_times (bool, optional) – If True, average spectra over the time axis. Default: False.

• fold (bool, optional) – Whether to fold the power spectrum in \(|k_\parallel|\). Default: False.

• plot_noise (bool, optional) – Whether to plot noise power spectrum curves or not. Default: False.

• delay (bool, optional) – Whether to plot the power spectrum in delay units (ns) or cosmological units (h/Mpc). Default: True.

• deltasq (bool, optional) – If True, plot dimensionless power spectra, Delta^2. This is ignored if delay=True. Default: False.

• legend (bool, optional) – Whether to switch on the plot legend. Default: False.

• ax (matplotlib.axes, optional) – Use this to pass in an existing Axes object, which the power spectra will be added to. (Warning: Labels and legends will not be altered in this case, even if the existing plot has completely different axis labels etc.) If None, a new Axes object will be created. Default: None.

• component (str, optional) – Component of complex spectra to plot, options=[‘abs’, ‘real’, ‘imag’]. Default: ‘real’.

• lines (bool, optional) – If True, plot lines between bandpowers for a given pspectrum. Default: True.

• markers (bool, optional) – If True, plot circles at bandpowers. Filled for positive, empty for negative. Default: False.

• error (str, optional) – If not None and if error exists in uvp stats_array, plot errors on bandpowers. Default: None.

• times (array_like, optional) – Float ndarray containing elements from time_avg_array to plot.

• logscale (bool, optional) – If True, put plot on a log-scale. Else linear scale. Default: True.

• force_plot (bool, optional) – If plotting a large number of spectra (>100), this function will error. Set this to True to override this large plot error and force plot. Default: False.

• label_type (int, optional) –

  Line label type in legend, options=[‘key’, ‘blpair’, ‘blpairt’].

  key : Label lines based on (spw, blpair, pol) key. blpair : Label lines based on blpair. blpairt : Label lines based on blpair-time.

• plot_stats (string, optional) – If not None, plot an entry in uvp.stats_array instead of power spectrum in uvp.data_array.

• kwargs (dict, optional) – Extra kwargs to pass to _all_ ax.plot calls.

Returns:

fig – Matplotlib Figure instance.

Return type:

matplotlib.pyplot.Figure

Indices and tables

• Index
• Module Index
• Search Page