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=None, wgts=None, dsets_std=None, labels=None, beam=None, cals=None, cal_flag=True)[source]
__init__(dsets=None, wgts=None, dsets_std=None, labels=None, beam=None, cals=None, cal_flag=True)[source]

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)[source]

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()[source]

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, exact_windows=False, ftbeam=None, verbose=True, filter_extensions=None, exact_norm=False, history='', r_params=None, cov_model='empirical', known_cov=None, allow_fft=False)[source]

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

  • exact_windows (bool, optional) – If True, compute exact window functions and sets store_window=True. Default: False

  • ftbeam (str or FTBeam, optional) – Definition of the beam Fourier transform to be used. Options include;

    • Root name of the file to use, without the polarisation

    Ex : FT_beam_HERA_dipole (+ path) - ‘’ for computation from beam simulations (slow) - FTBeam object. Make sure polarisations and bandwidths are consistent with the data set.

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

  • allow_fft (bool, optional) – Use an fft to compute q-hat. Default is False.

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)[source]

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)[source]

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)[source]

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[source]

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

__init__()[source]

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='')[source]

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_avgbool, optional

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

blpair_weightslist, 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.

inplacebool, optional

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

add_to_historystr, optional

Added text to add to file history.

Notes

Currently, every baseline-pair in a blpair group must have the same Ntpairs, 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()[source]

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', scalar=None)[source]

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=(Ntpairs,)

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

  • scalar (float) – Optional noise scalar used to convert from Tsys to cosmological units output from pspecbeam.compute_pspec_scalar(…,noise_scalar=True) Default is None. If None provided, will calculate scalar in function.

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)[source]

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()[source]

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=(Nbltpairs,)

Return type

float ndarray

get_data(key, omit_flags=False)[source]

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)[source]

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)[source]

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)[source]

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)[source]

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)[source]

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)[source]

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)[source]

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)[source]

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)[source]

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.