Support input data with the alternative baseline layout.

The MWA Rust libraries subscribe to the "antenna1 - antenna2" baseline layout,
but others like to use "antenna2 - antenna1" instead. This manifests in uvfits
and MS data as the negative UVW coordinates of what we would normally
generate. Not detecting and altering behaviour based off of this baseline layout
difference means that any simulated visibilities are reflected in the UV
plane.

When the alternative baseline order is detected, any visibilities that are read
are conjugated.

The baseline layout is detected by (1) finding the first UVW coordinates within
the data, (2) comparing it to the UVW we would generate from the provided
XYZs, (3) comparing it to the UVW we would generate from the XYZs if they were
negated, and (4) using the smallest difference in UVW to inform what the
baseline layout is.

CHANGELOG.md

@@ -8,6 +8,7 @@ Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [0.2.2] - Unreleased
 ### Added
+- Support for visibilities using the ant2-ant1 ordering rather than ant1-ant2.
 - Add new errors
   - If all baselines are flagged due to UVW cutoffs, then the program exits with
     this report.

src/io/read/mod.rs

@@ -18,7 +18,11 @@ pub use uvfits::UvfitsReader;
 
 use std::collections::HashSet;
 
-use marlu::{Jones, MwaObsContext as MarluMwaObsContext};
+use hifitime::{Duration, Epoch};
+use marlu::{
+    precession::precess_time, Jones, LatLngHeight, MwaObsContext as MarluMwaObsContext, RADec,
+    XyzGeodetic, UVW,
+};
 use mwalib::MetafitsContext;
 use ndarray::prelude::*;
 use vec1::Vec1;
@@ -105,3 +109,57 @@ pub struct AutoData<'a, 'b, 'c> {
     pub weights_fb: ArrayViewMut2<'b, f32>,
     pub tile_baseline_flags: &'c TileBaselineFlags,
 }
+
+/// With a dataset's UVW and the XYZs that correspond to it, compare with a UVW
+/// that we form from the XYZs. This allows us to determine the "baseline order"
+/// that the software that wrote this dataset used. `hyperdrive` and friends use
+/// ant1-ant2, but others may use ant2-ant1. If we detect ant2-ant1, we know
+/// that we have to conjugate this dataset's visibilities if want to continue
+/// using ant1-ant2.
+///
+/// It is not anticipated that precession has an impact here.
+fn baseline_convention_is_different(
+    data_uvw: UVW,
+    tile1_xyz: XyzGeodetic,
+    tile2_xyz: XyzGeodetic,
+    array_position: LatLngHeight,
+    phase_centre: RADec,
+    first_timestamp: Epoch,
+    dut1: Option<Duration>,
+) -> bool {
+    let precession_info = precess_time(
+        array_position.longitude_rad,
+        array_position.latitude_rad,
+        phase_centre,
+        first_timestamp,
+        dut1.unwrap_or_default(),
+    );
+    let xyzs = precession_info.precess_xyz(&[tile1_xyz, tile2_xyz]);
+    let UVW {
+        u: u_p1,
+        v: v_p1,
+        w: w_p1,
+    } = UVW::from_xyz(
+        xyzs[0] - xyzs[1],
+        phase_centre.to_hadec(precession_info.lmst_j2000),
+    );
+    let UVW {
+        u: u_p2,
+        v: v_p2,
+        w: w_p2,
+    } = UVW::from_xyz(
+        xyzs[1] - xyzs[0],
+        phase_centre.to_hadec(precession_info.lmst_j2000),
+    );
+
+    // Which UVW is closer to the data?
+    let UVW { u, v, w } = data_uvw;
+    let diff1 = (u - u_p1).abs() + (v - v_p1).abs() + (w - w_p1).abs();
+    let diff2 = (u - u_p2).abs() + (v - v_p2).abs() + (w - w_p2).abs();
+
+    // If `diff2` is smaller than `diff1`, then the standard baseline order is
+    // good, no need to do anything else. Otherwise, the other baseline order is
+    // used; the tile XYZs need to be negated and the visibility data need to be
+    // complex conjugated.
+    diff2 < diff1
+}

src/io/read/ms/mod.rs

@@ -115,6 +115,10 @@ pub struct MsReader {
     /// Is the weight column two-dimensional? We track this here because it
     /// could be 1D or 2D.
     weight_col_is_2d: bool,
+
+    /// If the incoming data uses ant2-ant1 UVWs instead of ant1-ant2 UVWs, we
+    /// need to conjugate the visibilities to match what will be modelled.
+    conjugate_vis: bool,
 }
 
 impl MsReader {
@@ -825,6 +829,52 @@ impl MsReader {
         }
         .map(Duration::from_seconds);
 
+        // Compare the first cross-correlation row's UVWs against UVWs that we
+        // would make with the existing tiles. If they're negative of one
+        // another, we need to negate our XYZs to match the UVWs the data use.
+        let conjugate_vis = {
+            let mut data_uvw = None;
+            let mut tile1_xyz = None;
+            let mut tile2_xyz = None;
+
+            for i_row in 0.. {
+                let uvw: Vec<f64> = main_table.get_cell_as_vec("UVW", i_row)?;
+                let ant1: i32 = main_table.get_cell("ANTENNA1", i_row)?;
+                let ant2: i32 = main_table.get_cell("ANTENNA2", i_row)?;
+
+                // We need to use a cross-correlation baseline.
+                if ant1 == ant2 {
+                    continue;
+                }
+
+                data_uvw = Some(UVW {
+                    u: uvw[0],
+                    v: uvw[1],
+                    w: uvw[2],
+                });
+                tile1_xyz = Some(tile_xyzs[ant1 as usize]);
+                tile2_xyz = Some(tile_xyzs[ant2 as usize]);
+                break;
+            }
+
+            if baseline_convention_is_different(
+                // If this data somehow has no cross-correlation data, using
+                // default values won't affect anything.
+                data_uvw.unwrap_or_default(),
+                tile1_xyz.unwrap_or_default(),
+                tile2_xyz.unwrap_or_default(),
+                array_position,
+                phase_centre,
+                *timestamps.first(),
+                dut1,
+            ) {
+                warn!("uvfits UVWs use the other baseline convention; will conjugate incoming visibilities");
+                true
+            } else {
+                false
+            }
+        };
+
         let obs_context = ObsContext {
             obsid,
             timestamps,
@@ -862,6 +912,7 @@ impl MsReader {
             data_col_name,
             weight_col_name,
             weight_col_is_2d,
+            conjugate_vis,
         };
         Ok(ms)
     }
@@ -1106,26 +1157,81 @@ impl MsReader {
 
         // Transform the data, depending on what the actual polarisations are.
         if let Some(crosses) = crosses.as_mut() {
-            match self.obs_context.polarisations {
-                // These are all handled correctly.
-                Polarisations::XX_XY_YX_YY => (),
-                Polarisations::XX => (),
+            let c0 = num_complex::Complex32::default();
+            match (self.conjugate_vis, self.obs_context.polarisations) {
+                // These pols are all handled correctly.
+                (false, Polarisations::XX_XY_YX_YY) => (),
+                (false, Polarisations::XX) => (),
+                // Just conjugate.
+                (true, Polarisations::XX_XY_YX_YY | Polarisations::XX) => {
+                    crosses.vis_fb.mapv_inplace(|j| {
+                        Jones::from([j[0].conj(), j[1].conj(), j[2].conj(), j[3].conj()])
+                    })
+                }
 
                 // Because we read in one polarisation, it was treated as XX,
                 // but this is actually YY.
-                Polarisations::YY => crosses.vis_fb.mapv_inplace(|j| {
-                    Jones::from([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, j[0].re, j[0].im])
-                }),
+                (false, Polarisations::YY) => crosses
+                    .vis_fb
+                    .mapv_inplace(|j| Jones::from([c0, c0, c0, j[0]])),
+                (true, Polarisations::YY) => crosses
+                    .vis_fb
+                    .mapv_inplace(|j| Jones::from([c0, c0, c0, j[0].conj()])),
+
                 // What looks like XY is YY.
-                Polarisations::XX_YY => crosses.vis_fb.mapv_inplace(|j| {
-                    Jones::from([j[0].re, j[0].im, 0.0, 0.0, 0.0, 0.0, j[1].re, j[1].im])
-                }),
+                (false, Polarisations::XX_YY) => crosses
+                    .vis_fb
+                    .mapv_inplace(|j| Jones::from([j[0], c0, c0, j[1]])),
+                (true, Polarisations::XX_YY) => crosses
+                    .vis_fb
+                    .mapv_inplace(|j| Jones::from([j[0].conj(), c0, c0, j[1].conj()])),
+
+                // What looks like YX is YY.
+                (false, Polarisations::XX_YY_XY) => crosses
+                    .vis_fb
+                    .mapv_inplace(|j| Jones::from([j[0], j[1], c0, j[2]])),
+                (true, Polarisations::XX_YY_XY) => crosses
+                    .vis_fb
+                    .mapv_inplace(|j| Jones::from([j[0].conj(), j[1].conj(), c0, j[2].conj()])),
+            }
+        }
+        if let Some(autos) = autos.as_mut() {
+            let c0 = num_complex::Complex32::default();
+            match (self.conjugate_vis, self.obs_context.polarisations) {
+                // These pols are all handled correctly.
+                (false, Polarisations::XX_XY_YX_YY) => (),
+                (false, Polarisations::XX) => (),
+                // Just conjugate.
+                (true, Polarisations::XX_XY_YX_YY | Polarisations::XX) => {
+                    autos.vis_fb.mapv_inplace(|j| {
+                        Jones::from([j[0].conj(), j[1].conj(), j[2].conj(), j[3].conj()])
+                    })
+                }
+
+                // Because we read in one polarisation, it was treated as XX,
+                // but this is actually YY.
+                (false, Polarisations::YY) => autos
+                    .vis_fb
+                    .mapv_inplace(|j| Jones::from([c0, c0, c0, j[0]])),
+                (true, Polarisations::YY) => autos
+                    .vis_fb
+                    .mapv_inplace(|j| Jones::from([c0, c0, c0, j[0].conj()])),
+
+                // What looks like XY is YY.
+                (false, Polarisations::XX_YY) => autos
+                    .vis_fb
+                    .mapv_inplace(|j| Jones::from([j[0], c0, c0, j[1]])),
+                (true, Polarisations::XX_YY) => autos
+                    .vis_fb
+                    .mapv_inplace(|j| Jones::from([j[0].conj(), c0, c0, j[1].conj()])),
+
                 // What looks like YX is YY.
-                Polarisations::XX_YY_XY => crosses.vis_fb.mapv_inplace(|j| {
-                    Jones::from([
-                        j[0].re, j[0].im, j[1].re, j[1].im, 0.0, 0.0, j[2].re, j[2].im,
-                    ])
-                }),
+                (false, Polarisations::XX_YY_XY) => autos
+                    .vis_fb
+                    .mapv_inplace(|j| Jones::from([j[0], j[1], c0, j[2]])),
+                (true, Polarisations::XX_YY_XY) => autos
+                    .vis_fb
+                    .mapv_inplace(|j| Jones::from([j[0].conj(), j[1].conj(), c0, j[2].conj()])),
             }
         }
 

src/io/read/ms/tests.rs

@@ -920,9 +920,9 @@ fn test_sdc3() {
 
     assert_abs_diff_eq!(
         cross_vis[(0, 0)],
-        // N.B. This should be conjugated, because the SDC3 data uses ant2-ant1
-        // UVWs, whereas we expect ant1-ant2.
-        Jones::from([-8.517027, 7.5777674, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
+        // N.B. This is conjugated from what's in the data, because the SDC3
+        // data uses ant2-ant1 UVWs, whereas we expect ant1-ant2.
+        Jones::from([-8.517027, -7.5777674, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
     );
     assert_abs_diff_eq!(cross_vis_weights[(0, 0)], 1.0);
 }