Merge pull request #37 from MWATelescope/fix_uvfits_baseline_decode

Fix uvfits baseline decode

Author: gsleap

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "marlu"
-version = "0.16.0"
+version = "0.16.1"
 authors = [
     "Christopher H. Jordan <christopherjordan87@gmail.com>",
     "Dev Null <dev.null@curtin.edu.au>",

RELEASES.md

@@ -1,5 +1,10 @@
 <!-- markdownlint-disable=MD025 -->
 
+# Version 0.16.1 (2025-04-11)
+
+- Fixed decode_uvfits_baseline bug (#36)
+- minor code reformatting
+
 # Version 0.16.0 (2025-04-10)
 
 - update mwalib 1.8.7

src/averaging.rs

@@ -203,7 +203,7 @@ pub type VisData33 = (Array3<Jones<f32>>, Array3<f32>);
 /// chunk of input visibilities:
 /// - unflagged weights are added together
 /// - if all visibilities in a chunk are flagged, then the result is the geometric
-///     mean of the chunk.
+///   mean of the chunk.
 /// - otherwise the visibility is the weighted mean of the unflagged visibilities.
 ///
 /// This has been validated thoroughly against Cotter.

src/io/ms.rs

@@ -464,7 +464,7 @@ impl MeasurementSetWriter {
     /// - `name` - Spectral Window name (`NAME` column)
     /// - `ref_freq` - Reference frequency (`REF_FREQUENCY` column)
     /// - `chan_info` - A two-dimensional array of shape (n, 4), containing the
-    ///     following for each channel:
+    ///   following for each channel:
     ///     - `CHAN_FREQ` - the center frequencies
     ///     - `CHAN_WIDTH` - channel widths,
     ///     - `EFFECTIVE_BW` - effective noise bandwidths
@@ -521,14 +521,14 @@ impl MeasurementSetWriter {
     /// - `name` - Spectral Window name (`NAME` column)
     /// - `ref_freq` - Reference frequency (`REF_FREQUENCY` column)
     /// - `chan_info` - A two-dimensional array of shape (n, 4), containing the
-    ///     following for each channel:
+    ///   following for each channel:
     ///     - `CHAN_FREQ` - the center frequencies
     ///     - `CHAN_WIDTH` - channel widths,
     ///     - `EFFECTIVE_BW` - effective noise bandwidths
     ///     - `RESOLUTION` - resolutions.
     /// - `total_bw` - Total bandwidth (`TOTAL_BANDWIDTH` column)
     /// - `centre_subband_nr` - This is the "sky" channel number of the center coarse channel in
-    ///     the spectral window.
+    ///   the spectral window.
     /// - `flag` - Row flag (`FLAG_ROW` column)
     #[allow(clippy::too_many_arguments)]
     pub fn write_spectral_window_row_mwa(
@@ -677,7 +677,7 @@ impl MeasurementSetWriter {
     /// - `idx` - row index to write to (ensure enough rows have been added)
     /// - `corr_type` - The polarization type for each correlation product, as a Stokes enum.
     /// - `corr_product` - Indices describing receptors of feed going into correlation.
-    ///     Shape should be [n, 2] where n is the length of `corr_type`
+    ///   Shape should be [n, 2] where n is the length of `corr_type`
     /// - `flag_row` - Row flag
     #[allow(clippy::ptr_arg)]
     pub fn write_polarization_row(
@@ -771,8 +771,8 @@ impl MeasurementSetWriter {
     /// - `code` - Special characteristics of field, e.g. Bandpass calibrator
     /// - `time` - Time origin for direction and rate
     /// - `dir_info` - An array of polynomial coefficients to calculate a direction
-    ///     (RA, DEC) relative to `time`. The shape is  [3, p, 2], where p is the maximum
-    ///     order of all polynomials, and there are three direction polynomials:
+    ///   (RA, DEC) relative to `time`. The shape is  [3, p, 2], where p is the maximum
+    ///   order of all polynomials, and there are three direction polynomials:
     ///     - `DELAY_DIR` - Direction of delay center (e.g. RA, DEC) in time
     ///     - `PHASE_DIR` - Direction of phase center (e.g. RA, DEC) in time
     ///     - `REFERENCE_DIR` - Direction of reference center (e.g. RA, DEC) in time
@@ -830,8 +830,8 @@ impl MeasurementSetWriter {
     /// - `code` - Special characteristics of field, e.g. Bandpass calibrator
     /// - `time` - Time origin for direction and rate
     /// - `dir_info` - An array of polynomial coefficients to calculate a direction
-    ///     (RA, DEC) relative to `time`. The shape is  [3, p, 2], where p is the maximum
-    ///     order of all polynomials, and there are three direction polynomials:
+    ///   (RA, DEC) relative to `time`. The shape is  [3, p, 2], where p is the maximum
+    ///   order of all polynomials, and there are three direction polynomials:
     ///     - `DELAY_DIR` - Direction of delay center (e.g. RA, DEC) in time
     ///     - `PHASE_DIR` - Direction of phase center (e.g. RA, DEC) in time
     ///     - `REFERENCE_DIR` - Direction of reference center (e.g. RA, DEC) in time
@@ -1605,7 +1605,7 @@ impl MeasurementSetWriter {
     /// - `state_id` - ID for this observing state
     /// - `sigma` - Estimated rms noise for channel with unity bandpass response
     /// - `data` - an `[n, p]` shaped ndarray of complex visibilities, where `n`
-    ///     is the number of channels, and p is the number of polarizations
+    ///   is the number of channels, and p is the number of polarizations
     /// - `flags` - an `[n, p]` shaped ndarray of boolean flags.
     /// - `weights` - a `[p]` shaped ndarray of weights for each polarization
     ///

src/io/uvfits.rs

@@ -80,7 +80,7 @@ pub const fn encode_uvfits_baseline(ant1: usize, ant2: usize) -> usize {
 /// start at 1.
 #[allow(dead_code)]
 pub const fn decode_uvfits_baseline(bl: usize) -> (usize, usize) {
-    if bl < 65_535 {
+    if bl <= 65_535 {
         let ant2 = bl % 256;
         let ant1 = (bl - ant2) / 256;
         (ant1, ant2)
@@ -2830,4 +2830,110 @@ mod tests {
             }
         }
     }
+
+    #[test]
+    fn test_encode_uvfits_baseline_256t() {
+        // test hand calculated edge cases
+        assert_eq!(encode_uvfits_baseline(1,1), 257);
+        assert_eq!(encode_uvfits_baseline(1,2), 258);
+        assert_eq!(encode_uvfits_baseline(1,255), 511);
+        assert_eq!(encode_uvfits_baseline(1,256), 67840);
+        assert_eq!(encode_uvfits_baseline(127,127), 32639);
+        assert_eq!(encode_uvfits_baseline(127,255), 32767);
+        assert_eq!(encode_uvfits_baseline(127,256), 325888);
+        assert_eq!(encode_uvfits_baseline(128,128), 32896);
+        assert_eq!(encode_uvfits_baseline(128,255), 33023);
+        assert_eq!(encode_uvfits_baseline(128,256), 327936);
+        assert_eq!(encode_uvfits_baseline(254,254), 65278);
+        assert_eq!(encode_uvfits_baseline(254,255), 65279);
+        assert_eq!(encode_uvfits_baseline(254,256), 585984);
+        assert_eq!(encode_uvfits_baseline(255,255), 65535);
+        assert_eq!(encode_uvfits_baseline(255,256), 588032);
+        assert_eq!(encode_uvfits_baseline(256,256), 590080);
+    }
+
+    #[test]
+    fn test_decode_uvfits_baseline_256t() {
+        // Test hand calculated edge cases
+        assert_eq!(decode_uvfits_baseline(257), (1,1));
+        assert_eq!(decode_uvfits_baseline(258), (1,2));
+        assert_eq!(decode_uvfits_baseline(511), (1,255));
+        assert_eq!(decode_uvfits_baseline(67840), (1,256));
+        assert_eq!(decode_uvfits_baseline(32639), (127,127));
+        assert_eq!(decode_uvfits_baseline(32767), (127,255));
+        assert_eq!(decode_uvfits_baseline(325888), (127,256));
+        assert_eq!(decode_uvfits_baseline(32896), (128,128));
+        assert_eq!(decode_uvfits_baseline(33023), (128,255));
+        assert_eq!(decode_uvfits_baseline(327936), (128,256));
+        assert_eq!(decode_uvfits_baseline(65278), (254,254));
+        assert_eq!(decode_uvfits_baseline(65279), (254,255));
+        assert_eq!(decode_uvfits_baseline(585984), (254,256));
+        assert_eq!(decode_uvfits_baseline(65535), (255,255));
+        assert_eq!(decode_uvfits_baseline(588032), (255,256));
+        assert_eq!(decode_uvfits_baseline(590080), (256,256));
+    }
 }
+
+#[test]
+fn test_encode_decode_uvfits_all_baselines_128t() {
+    // Test all baselines can be encoded then decoded correctly    
+    // MWAX 128T
+    let n_ants: usize= 128;
+    
+    for ant1 in 1..=n_ants {
+        for ant2 in ant1..=n_ants {
+            let bl = encode_uvfits_baseline(ant1, ant2);
+            let (a1,a2) = decode_uvfits_baseline(bl);
+            assert_eq!(a1, ant1);
+            assert_eq!(a2, ant2);
+        }
+    }
+}
+
+#[test]
+fn test_encode_decode_uvfits_all_baselines_256t() {
+    // Test all baselines can be encoded then decoded correctly
+    // MWAX 256T
+    let n_ants: usize= 256;
+    
+    for ant1 in 1..=n_ants {
+        for ant2 in ant1..=n_ants {
+            let bl = encode_uvfits_baseline(ant1, ant2);
+            let (a1,a2) = decode_uvfits_baseline(bl);
+            assert_eq!(a1, ant1);
+            assert_eq!(a2, ant2);
+        }
+    }
+}
+
+#[test]
+fn test_encode_decode_uvfits_all_baselines_264t() {
+    // Test all baselines can be encoded then decoded correctly
+    // MWAX 256T + 1 more receiver
+    let n_ants: usize= 264;
+    
+    for ant1 in 1..=n_ants {
+        for ant2 in ant1..=n_ants {
+            let bl = encode_uvfits_baseline(ant1, ant2);
+            let (a1,a2) = decode_uvfits_baseline(bl);
+            assert_eq!(a1, ant1);
+            assert_eq!(a2, ant2);
+        }
+    }
+}
+
+#[test]
+fn test_encode_decode_uvfits_all_baselines_512t() {
+    // Test all baselines can be encoded then decoded correctly
+    // MWAX 512T (Phase 4(!))
+    let n_ants: usize= 512;
+    
+    for ant1 in 1..=n_ants {
+        for ant2 in ant1..=n_ants {
+            let bl = encode_uvfits_baseline(ant1, ant2);
+            let (a1,a2) = decode_uvfits_baseline(bl);
+            assert_eq!(a1, ant1);
+            assert_eq!(a2, ant2);
+        }
+    }
+}