Birli

A Rust library for common preprocessing tasks performed in the data pipeline of the Murchison Widefield Array (MWA), located on the land of the Wajarri Yamatji people in Murchison Shire, Western Australia.

Birli reads MWA correlator visibilities in the gpufits file format using mwalib, which supports the existing "legacy" MWA correlator, as well as the in-development "MWAX" correlator.

Birli is the Wajarri word for lightning, a common cause of outages at the MWA, and a great descriptor for the speed which this library intends to deliver.

Installation

Prerequisites

• A Rust compiler with a version >= 1.65.0 - https://www.rust-lang.org/tools/install
• AOFlagger >= 3.0 (Ubuntu > 21.04: apt install aoflagger-dev)
• CFitsIO >= 3.49 (Ubuntu > 20.10: apt install libcfitsio-dev)

for OS-specific instructions, check out the linux CI Script; the Makefile.toml; and the Dockerfile as these are tested regularly. The instructions below may be updated less frequently, but are better documented.

(Debian/Ubuntu) Linux Setup

# Prerequisites for rustup, cargo and cargo-make
sudo apt install -y gcc libssl-dev pkg-config curl unzip wget
# Run the Rustup install script, profile=default, toolchain=stable
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs -sSf | sh -s -- -y
# Cargo make uses Makefile.toml to automate development tasks
cargo install --force cargo-make
# Use multiple cores when compiling C/C++ libraries
export MAKEFLAGS="-j $MAKEFLAGS"
# Install prerequisite C/C++ libraries
cargo make install_deps
# Ensure that rust can find the C/C++ libraries.
# AOFlagger and CFitsIO default to /usr/local/lib
export LD_LIBRARY_PATH="/usr/local/lib/"

macOS

Previously macOS was supported, however this has ben dropped due stability issues with casacore and aoflagger.

The following may or may not work.

brew tap mwatelescope/tap
brew install mwatelescope/tap/aoflagger
AOFLAGGER_INCLUDE_DIR=/opt/homebrew/include/ cargo install --path .
export DYLD_FALLBACK_LIBRARY_PATH=/opt/homebrew/lib/
birli

It is not possible to set the library search path globally https://github.com/Homebrew/brew/issues/13481. So if you don't set DYLD_FALLBACK_LIBRARY_PATH you will see an error like

libc++abi: terminating due to uncaught exception of type std::runtime_error: cannot open : No such file or directory

in some cases, macOS will deliberately unset this.

export DYLD_FALLBACK_LIBRARY_PATH=/opt/homebrew/lib/
echo ${DYLD_FALLBACK_LIBRARY_PATH:-unset} # it's set here
bash -c 'echo ${DYLD_FALLBACK_LIBRARY_PATH:-unset}' # but not here

Other Operating Systems

Unfortunately most of the prerequisites aren't available on Windows. However, WSL is great, and there is a docker image! You could use VSCode remote for WSL or Docker. Your best best is Ubuntu LTS

Installing the binary

cargo install --path .

This creates a birli binary in $HOME/.cargo/bin

Troubleshooting

Test suite

Having issues with Birli? run the test suite to narrow down your issue.

cargo test

Dependencies

Experiencing segfaults? I can guarantee it's because of one of the C library dependencies. Make sure you have the right versions of all the libraries. These are specified in Prerequisites.

Get library versions on linux with:

pkg-config --modversion cfitsio
aoflagger --version

If you have something like CASA installed from apt, it's going to put an ancient cfitsio library version in /usr/lib/x86_64-linux-gnu/, to get around this, you must export LD_LIBRARY_PATH=/usr/local/lib/ in the shell so that Birli can find the correct library version.

Logging

You can enable additional logging on individual Rust modules by setting the RUST_LOG environment variable. For example:

RUST_LOG=trace birli ... # set log level to trace for all module (including dependencies)
RUST_LOG=birli=debug birli ... # set log level to debug for birli only
RUST_LOG=birli::io=error birli ... # only show warnings for birli's io module

For more examples, see the env_logger docs

The default log level in info

Docker

Couldn't get it working on your environment? You can always run Birli in Docker

docker run mwatelescope/birli:latest -h

Want to open a shell within a fully provisioned Birli development environment? Easy!

docker run -it --entrypoint /bin/bash --volume $PWD:/app mwatelescope/birli:latest

Note: This mounts the current directory to /app in the Docker image, meaning both of these systems share the same target folder. so if your host system is a different architecture than Docker, you may need to cargo clean each time you switch between these environments. You may also want to temporarily disable any linters or language servers that use

Singularity on HPC

# - load the singularity module
module load singularity
# - cd into your preferred sif file location, e.g. /pawsey/mwa/singularity/birli
# - create a .sif file from the latest mwatelescope/birli docker image
singularity pull --dir . docker://mwatelescope/birli:latest
# - run birli within the singularity image
singularity exec  /pawsey/mwa/singularity/birli/birli_latest.sif /app/target/release/birli ${YOUR_BIRLI_ARGS}

see this gist for an example of a Garrawarla SLURM job using Birli

Singularity on HPC (debug mode)

This will give you much more information about any problem you're having with Birli, however the debug build is not optimised, and is much slower.

# - request an interactive HPC session
salloc --partition workq --time 1:00:00 --nodes 1 -c 38 --mem=350G
# - load the singularity module
module load singularity
# - cd into your preferred sif file location, e.g. /pawsey/mwa/singularity/birli
# - create a .sif file from the latest mwatelescope/birli docker image
singularity pull --dir . docker://mwatelescope/birli:debug
# - run birli within the singularity image
singularity exec  /pawsey/mwa/singularity/birli/birli_debug.sif /bin/bash

then within this shell

# - enable lots of logs
export RUST_LOG=trace
# - run birli in debug mode with GDB
gdb --args /app/target/debug/birli ${YOUR_BIRLI_ARGS}
# > run

Usage

birli -h

USAGE:
    birli [OPTIONS] --metafits <PATH> <PATHS>...

OPTIONS:
        --apply-di-cal <PATH>        Apply DI calibration solutions before averaging
        --dry-run                    Just print the summary and exit
        --emulate-cotter             Use Cotter's array position, not MWAlib's
    -h, --help                       Print help information
        --no-draw-progress           do not show progress bars
        --phase-centre <RA> <DEC>    Override Phase centre from metafits (degrees)
        --pointing-centre            Use pointing instead phase centre
    -V, --version                    Print version information

INPUT:
    -m, --metafits <PATH>    Metadata file for the observation
    <PATHS>...           GPUBox files to process

SELECTION:
        --no-sel-autos                [WIP] Deselect autocorrelations
        --no-sel-flagged-ants         [WIP] Deselect flagged antennas
        --provided-chan-ranges        Only consider provided channels
        --sel-ants <ANTS>...          [WIP] Antenna to select
        --sel-chan-ranges <RANGES>    Select separate channel ranges
        --sel-time <MIN> <MAX>        Timestep index range (inclusive) to select

RESOURCE LIMITS:
        --max-memory <GIBIBYTES>    Estimate --time-chunk with <GIBIBYTES> GiB each chunk.
        --time-chunk <STEPS>        Process observation in chunks of <STEPS> timesteps.

FLAGGING:
        --flag-antennas <ANTS>...         Flag antenna indices
        --flag-autos                      Flag auto correlations
        --flag-coarse-chans <CHANS>...    Flag additional coarse chan indices
        --flag-dc                         Force flagging of DC centre chans
        --flag-edge-chans <COUNT>         Flag <COUNT> fine chans on the ends of each coarse
        --flag-edge-width <KHZ>           Flag bandwidth [kHz] at the ends of each coarse chan
        --flag-end <SECONDS>              Flag seconds before the last provided time
        --flag-end-steps <COUNT>          Flag <COUNT> steps before the last provided
        --flag-fine-chans <CHANS>...      Flag fine chan indices in each coarse chan
        --flag-init <SECONDS>             Flag <SECONDS> after first common time (quack time)
        --flag-init-steps <COUNT>         Flag <COUNT> steps after first common time
        --flag-times <STEPS>...           Flag additional time steps
        --no-flag-dc                      Do not flag DC centre chans
        --no-flag-metafits                Ignore antenna flags in metafits

CORRECTION:
        --no-cable-delay           Do not perform cable length corrections
        --no-digital-gains         Do not perform digital gains corrections
        --no-geometric-delay       Do not perform geometric corrections
        --passband-gains <TYPE>    Type of PFB passband filter gains correction to apply [default:
                                   auto] [possible values: none, cotter, jake, auto]
        --van-vleck                Apply Van Vleck corrections

AVERAGING:
        --avg-freq-factor <FACTOR>    Average <FACTOR> channels per averaged channel
        --avg-freq-res <KHZ>          Frequency resolution of averaged data
        --avg-time-factor <FACTOR>    Average <FACTOR> timesteps per averaged timestep
        --avg-time-res <SECONDS>      Time resolution of averaged data

OUTPUT:
    -f, --flag-template <TEMPLATE>    The template used to name flag files. Percents are substituted
                                      for the zero-prefixed GPUBox ID, which can be up to 3
                                      characters long. Example: FlagFile%%%.mwaf
    -M, --ms-out <PATH>               Path for measurement set output
    -u, --uvfits-out <PATH>           Path for uvfits output

AOFLAGGER:
        --aoflagger-strategy <PATH>    Strategy to use for RFI Flagging
        --no-rfi                       Do not perform RFI Flagging with aoflagger

Note: the aoflagger options are only available when the aoflagger feature is enabled.

Correction Details

Operations are performed in the order described by the following sections.

flowchart TD;
classDef in fill:#2aa198;
classDef out fill:#d33682;
classDef file fill:#268bd2;
classDef proc fill:#b58900;
classDef decision fill:#cb4b16;

metafits[/ metafits /]; class metafits file;
raw[/ raw data /]; class raw file;
prepUVFits[/ uvfits /]; class prepUVFits file;
prepMS[/ CASA Measurement Set /]; class prepMS file;
calSols[/"Calibration Solutions (.bin)"/]; class calSols file;
mwaf[/"flags"/]; class mwaf file;

vv[["Van Vleck
(iff <code>--van-vleck</code>)"]]; class vv proc;
cableDel[["Cable Delays (unless
<code>--no-cable-delay</code>)"]]; class cableDel proc;
digGains[["Digital Gains (unless
<code>--no-digital-gains</code>)"]]; class digGains proc;
PFB[["PFB Passband
(<code>--passband-gains</code>)"]]; class PFB proc;
AOFlag[["AOFlagger (see
<code>--aoflagger-strategy</code>
and <code>--no-rfi</code>)"]]; class AOFlag proc;
applyCal[[ Apply Solutions ]]; class applyCal proc;
geometric[["Geometric Delays (unless
<code>--no-geometric-delay</code>)"]]; class geometric proc;
avg[["Averaging (see <code>--avg-*</code>)"]]; class avg proc;

metafits -->|"--metafits (-m)"| vv
raw --> vv
vv --> cableDel --> digGains --> PFB --> AOFlag ---> geometric
AOFlag -->|"--flag-template (-f)"| mwaf
calSols -->|"--apply-di-cal"| applyCal
geometric --> applyCal --> avg
avg -->|"--uvfits-out (-u)"| prepUVFits
avg -->|"--ms-out (-M)"| prepMS

Van Vleck Corrections

Van Vleck corrections are applied to the visibilities to correct for the effects of quantization in the legacy Correlator, resulting in a smoother bandpass, and a reduced false positive flag rate in AOFlagger.

These corrections are computationally expensive, and not needed for most science cases. They are most applicable to Epoch of Reionization science, where a flat bandpass is desired.

Birli's implementation borrows heavily from Pyxie Star's pyuvdata implementation described in Van Vleck Memo A, which is based on the formulation in Benkevitch 2016.

The current implementation only accounts for the 4+4 bit quantization of the legacy correlator, and does not account for the 5+5 bit quantization of the RRI receiver.

Birli solves the full integral with maximum precision, and does not use the Chebychev approximation described in the memo.

Cable Delay Corrections

Cable delay correction involves adjusting visibility phases to correct for the differences in electrical length of the cable between each tile and it's receiver.

Legacy MWA correlator observations do not typically have cable delays applied, however MWAX observations can. The CABLEDEL key in the metafits describes what geometric delays have been applied.

By default, Birli will apply cable length corrections. You can use --no-cable-delay to disable this.

A baseline's cable lengths are determined by the difference between a baseline's rfInput electrical lengths, as specified the the TILEDATA HDU of the metafits. Complex visibilities are phase-shifted by an angle determined by the electrical length, and the channel's frequency.

let angle = -2.0 * PI * electrical_length_m * freq_hz / SPEED_OF_LIGHT_IN_VACUUM_M_PER_S;

Digital Gain Corrections

Each input in the raw data is scaled by a factor for each coarse channel. This is defined in the metafits primary hdu in the Gains column. Birli corrects these digital gains by default, you can disable this with --no-digital-gains

Coarse PFB Passband Corrections

There are four options for correcting for the coarse channel polyphase filter bank (PFB) passband filter gains. cotter is best for legacy correlator observations, and jake is best for MWAX correlator observations. auto picks the most appropriate option based on information from the metafits, and none disables this correction.

Although the PFB passband filter shape is a function of the filter coefficients used in the receiver, it is also influenced by the fine channelization technique used, which is slightly different between the two correlators.

The MWAX channeliser is just an FFT with a rectangular window for its windowing function, of which, the ultra fine channels are then aggregated to your frequency resolution; while the legacy correlator had a 2nd stage fine-PFB with minimum bin width of 10kHz. The two different methods have a different amount of leakage between fine channels, which can effect the band shape across a coarse channel. More leakage leads to the band shape being smoothed out in frequency.

The cotter gains were source from Cotter's _sb128ChannelSubbandValue2014FromMemo in subbandpassband.cpp. These sometimes referred to as Levine Gains. Since these gains were computed at the base legacy correlator resolution of 10KHz, they will not work on resolutions that are not a multiple of 10kHz.

The jake gains (credit to Jake Jones) are described in this wiki article

When applying pfb gains to an observation that is not at the same resolution as the gains, the gains need to be averaged to fit the data, and the exact details of this averaging depends on the correlator type. For more dtails, see the mwa wiki on averaging fine channels

RFI Flagging

By default, Birli will flag the data using the default MWA strategy in AOFlagger. You can use the --no-rfi option to disable this, or the --aoflagger-strategy option to proived your own strategy file.

Geometric Delay Corrections (AKA Phase Tracking)

Geometric correction involves adjusting visibility phases to correct for the differences in distance that light from the phase center has to travel to reach each tile.

Legacy MWA correlator observations are not typically phase tracked, however MWAX observations can have phase tracking applied. The GEODEL card in the metafits describes what geometric delays have been applied.

By default, Birli will apply geometric corrections at the phase center if they have not already been applied. It determines the observations phase center from the RAPHASE and DECPHASE cards in the metafits. If these are not available, the pointing center cards (RA and DEC) from the metafits are used. You can use --no-geometric-delay to disable geometric corrections, as well as the --phase-centre and --pointing-centre options to override the phase center.

A baseline's geometric length is determined by the w component of it's UVW fourier-space vector, after applying precession and nutation to it's tiles' positions and the phase center to the J2000 epoch, accounting for stellar aberration. Complex visibilities are phase-shifted by an angle determined by the w-component, and the channel's frequency.

let angle = -2.0 * PI * uvw.w * freq_hz / SPEED_OF_LIGHT_IN_VACUUM_M_PER_S;

Calibration

Birli can apply direction independent calibration solutions using the --apply-di-cal flag. Solutions are applied before averaging. The number of channels in the un-averaged visibilities must be an integer multiple of the number of channels in the calibration solutions file. Unlike Cotter, Birli will handle calibration solutions where a NaN value is present by flagging any visibilities where a NaN is present.

Currently, only the MWA aocal format (.bin), historically generated by the calibrate binary in the mwa-reduce package is supported. This format is described here, however due to the ambiguous definition of the startTime and endTime fields, their values are ignored and so only a single timeblock of solutions can be applied.

Cotter Emulation

The --emulate-cotter flag ensures that outputs match Cotter as much as possible. You should only use this flag if you need to perform a direct comparison with Cotter.

By default, Birli will use the MWA array position from MWALib in order to calculate UVWs and geometric corrections. This is more accurate than the one that Cotter uses, and is the main source of error when doing direct comparisons.

This flag is used as part of the tests in src/main.rs to validate that Birli's output matches that of Cotter to within an acceptable margin.

Averaging

To average the data in time or frequency by a given whole number factor, you can provide the --avg-time-factor or --avg-freq-factor options. This can also be achieved with the --avg-time-res and --avg-freq-res options which take a duration [seconds] or ammount of bandwidth [kHz] respectively. This second group of options will choose the closest whole number averaging factor based on the resolution of the input data.

Output

Birli can output visibility data to uvfits or measurement set with --ms-out (-M) or --uvfits-out (-u). It can also output flags for each coarse channel in .mwaf format with --flag-template (-f), where the % characters in the template argument are replaced with the same zero-prefixed coarse channel identifiers that are used to identify the coarse channel GPUBox files that the coarse channel data came from. For legacy data, use two percentage characters, since the coarse channel identifier is the GPUBox number. However, for MWAX data, the coarse channel identifier is the channel number, which needs three digits.

Example: automatically determine flag template

export flag_template='Flagfile%%.mwaf'
if [ ${obsid} -gt 1300000000 ]; then
    flag_template='Flagfile_ch%%%.mwaf'
fi
birli \
  -f $flag_template \
  ...

When processing a set of coarse channels which are not contiguous in receiver channel number, a suffix will be added to the measurement set and uvfits filenames which indicates the coarse channel, or coarse channel range in that file.

Comparison with Cotter

The following table shows how Birli options map onto Cotter options:

Birli performs all the same default preprocessing steps as Cotter when no flags are provided. The exceptions are that we have not yet implemented flagging of auto-correlations, pruning of flagged antennas. This means that birli <in/out args> is equivalent to:

 cotter \
  -allowmissing \
  -noantennapruning \
  -noflagautos \
  -nostats \
  -flag-strategy <mwa default aoflagger strategy>
  <in/out args>

There is no intention of replicating the following options Birli at this point, so please open an issue if these are important to you:

• Coarse channel selection (-sbcount, -sbstart): This can be done by simply changing which coarse channel files are given in the CLI arguments)
• Dysco compression (-use-dysco, -dysco-config)
• Manual metadata specification (-a, -h, -i): This information is readily available from metafits.
• -offline-gpubox-format
• Quality statistics (-saveqs, -histograms, -skipwrite, -nostats)
• -noflagmissings: If an HDU is missing, it should always be flagged.
• -apply: only -full-apply is supported.
• -noalign: gpuboxes are always aligned.
• CPU limit (-j): Birli uses crossbeam for concurrency which intelligency uses the compute resources available. Strict resource limits can be achieved with cgroups.
• Memory percentage limit (-mem): Only -absmem is supported. Determining memory limits on HPC systems is unreliable, so we recommend manually specifying a memory limit instead.
• -sbpassband <file>
• -flagfiles <name> apply existing flags

Example: RFI Flagging, corrections, averaging, output

In this example, we use the aoflagger subcommand to:

• Perform RFI flagging using the MWA-default flagging strategy
• Perform geometric and cable length corrections
• average the data to 4 seconds, 160khz
• Output visibilities to .uvfits (-u)

birli \
  -m tests/data/1254670392_avg/1254670392.metafits \
  -f "/tmp/Flagfile.Birli.MWA.%%.mwaf" \
  -u "/tmp/1254670392.birli.uvfits" \
  --avg-time-res 4 --avg-freq-res 160 \
  tests/data/1254670392_avg/1254670392_*gpubox*.fits

The equivalent Cotter commands would be:

# output uvfits
cotter \
  -m tests/data/1254670392_avg/1254670392.metafits \
  -o "tests/data/1254670392_avg/1254670392.cotter.uvfits" \
  -allowmissing \
  -noantennapruning \
  -noflagautos \
  -nostats \
  -timeres 4 \
  -freqres 160 \
  -flag-strategy /usr/local/share/aoflagger/strategies/mwa-default.lua \
  tests/data/1254670392_avg/1254670392_20191009153257_gpubox*.fits

Contributing

Pull requests are welcome! Please do your best to ensure that the high standards of test coverage are maintained.

Before each commit, use cargo make ci to ensure your code is formatted correctly.

Acknowledgement

This scientific work uses data obtained from the Murchison Radio-astronomy Observatory. We acknowledge the Wajarri Yamatji people as the traditional owners of the Observatory site.

Coverage

This repo is approved by...