Merge pull request #62 from MWATelescope/mwax_updates

Mwax updates: Fixes #13, #59, #60, #61

.github/workflows/coverage.yml

@@ -11,7 +11,7 @@ jobs:
 
     steps:
     - name: Checkout sources
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
 
     - name: Install dependencies
       run: |

.github/workflows/releases.yml

@@ -16,7 +16,7 @@ jobs:
 
     steps:
     - name: Checkout sources
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
 
     - name: Install stable toolchain
       uses: actions-rs/toolchain@v1

.github/workflows/test_linux.yml

@@ -14,7 +14,7 @@ jobs:
 
     steps:
     - name: Checkout sources
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
 
     - name: Install stable toolchain
       uses: actions-rs/toolchain@v1

.github/workflows/test_macos.yml

@@ -14,7 +14,7 @@ jobs:
 
     steps:
     - name: Checkout sources
-      uses: actions/checkout@v2
+      uses: actions/checkout@v3
 
     - name: Install stable toolchain
       uses: actions-rs/toolchain@v1

.gitignore

@@ -14,3 +14,6 @@ include/
 /test_files/1101503312_1_timestep/*.sub
 /test_files/1101503312_1_timestep/*.dat
 ccov.zip
+/test_files/1101503312_mwax_vcs/*.sub
+/test_files/1101503312_vcs/*.dat
+/test_files/1370755832_mwax_vcs_os/*.sub

CHANGELOG.md

@@ -2,6 +2,14 @@
 
 Changes in each release are listed below.
 
+## 0.17.0 21-Aug-2023
+
+* Bumped MSRV to 1.64.
+* MetafitsContext now supports: OVERSAMP from metafits. Exposed as a bool: `oversampled`. If this is true, then the observation was taken using oversampled coarse channels.
+* VoltageContext now supports handling reading of oversampled subfiles.
+* MetafitsContext now supports: DR_FLAG and DR_PARAM from metafits. `deripple_applied`is a boolean and `dreipple_param` is a String. If deripple is true then a deripple has been applied to each coarse channel to smooth the passband.
+* Weights (packet occupancy) is now available via the CorrelatorContext struct method `read_weights_by_baseline()` and `read_weights_by_baseline_into_buffer()`.
+
 ## 0.16.4 16-Jun-2023
 
 * Modified rf_input.vcs_order to return `input` if `input` is > 255
@@ -25,10 +33,10 @@ Changes in each release are listed below.
 * FFI/C: `MetafitsMetadata`->`dec_phase_center_deg` is now NaN not 0 when missing from metafits file.
 * `MetafitsContext::grid_name` defaults to 'NOGRID' when key is missing from metafits file.
 * `MetafitsContext::grid_number` defaults to 0 when key is missing from metafits file.
-* `MetafitsContext::sun_alt_deg` is now an Option<f64> to handle cases where key is missing from metafits file. Via FFI/C it is represented by a NaN.
-* `MetafitsContext::sun_distance_deg` is now an Option<f64> to handle cases where key is missing from metafits file. Via FFI/C it is represented by a NaN.
-* `MetafitsContext::moon_distance_deg` is now an Option<f64> to handle cases where key is missing from metafits file. Via FFI/C it is represented by a NaN.
-* `MetafitsContext::jupiter_distance_deg` is now an Option<f64> to handle cases where key is missing from metafits file. Via FFI/C it is represented by a NaN.
+* `MetafitsContext::sun_alt_deg` is now an Option&lt;f64&gt; to handle cases where key is missing from metafits file. Via FFI/C it is represented by a NaN.
+* `MetafitsContext::sun_distance_deg` is now an Option&lt;f64&gt; to handle cases where key is missing from metafits file. Via FFI/C it is represented by a NaN.
+* `MetafitsContext::moon_distance_deg` is now an Option&lt;f64&gt; to handle cases where key is missing from metafits file. Via FFI/C it is represented by a NaN.
+* `MetafitsContext::jupiter_distance_deg` is now an Option&lt;f64&gt; to handle cases where key is missing from metafits file. Via FFI/C it is represented by a NaN.
 * Update module docs.
 * Update dependencies.
 * Fixed a bunch of clippy lints.
@@ -118,7 +126,7 @@ Changes in each release are listed below.
 
 ## 0.9.0 09-Aug-2021
 
-* Added mwa_version <Option<MWAVersion>> to MetafitsContext struct.
+* Added mwa_version &lt;Option&lt;MWAVersion&gt;&gt; to MetafitsContext struct.
 * When working only with a MetafitsContext, a None can be passed in lieu of an MWAVersion, and mwalib will attempt to determine the correct MWAVersion based on the MODE keyword from the metafits file.
 * Added method get_expected_volt_filename() function to MetafitsContext.
 * Added digital_gains, dipole_gains and dipole_delays to FFI rfinput struct.

Cargo.toml

@@ -1,13 +1,13 @@
 [package]
 name = "mwalib"
-version = "0.16.4"
+version = "0.17.0"
 homepage = "https://github.com/MWATelescope/mwalib"
 repository = "https://github.com/MWATelescope/mwalib"
 readme = "README.md"
 authors = ["Greg Sleap <greg.sleap@curtin.edu.au>",
            "Christopher H. Jordan <christopherjordan87@gmail.com>"]
 edition = "2021"
-rust-version = "1.60"
+rust-version = "1.64"
 description = "A library to simplify reading Murchison Widefield Array (MWA) raw visibilities, voltages and metadata."
 license = "MPL-2.0"
 keywords = ["radioastronomy", "mwa", "astronomy"]
@@ -43,13 +43,13 @@ clap = { version = "3.0.0", features = ["derive"], optional = true }
 env_logger = { version = "0.9.0", optional = true }
 
 [dev-dependencies]
-csv = "1.1.0"
-float-cmp = "0.9.0"
-tempdir = "0.3.6"
+csv = "1.2.*"
+float-cmp = "0.9.*"
+tempdir = "0.3.*"
 
 [build-dependencies]
-built = "0.5.0"
-cbindgen = { version = "0.24.0", default_features = false }
+built = "0.6.*"
+cbindgen = { version = "0.24.*", default_features = false }
 
 [[example]]
 name = "mwalib-data-dump"

src/correlator_context/mod.rs

@@ -94,8 +94,10 @@ pub struct CorrelatorContext {
 
     /// The number of bytes taken up by a scan/timestep in each gpubox file.
     pub num_timestep_coarse_chan_bytes: usize,
-    /// The number of floats in each gpubox HDU.
+    /// The number of floats in each gpubox visibility HDU.
     pub num_timestep_coarse_chan_floats: usize,
+    /// The number of floats in each gpubox weights HDU.
+    pub num_timestep_coarse_chan_weight_floats: usize,
     /// This is the number of gpubox files *per batch*.
     pub num_gpubox_files: usize,
     /// `gpubox_batches` *must* be sorted appropriately. See
@@ -338,6 +340,9 @@ impl CorrelatorContext {
             _ => Vec::new(),
         };
 
+        let weight_floats: usize =
+            metafits_context.num_baselines * metafits_context.num_visibility_pols;
+
         Ok(CorrelatorContext {
             metafits_context,
             mwa_version: gpubox_info.mwa_version,
@@ -373,6 +378,7 @@ impl CorrelatorContext {
             gpubox_time_map: gpubox_info.time_map,
             num_timestep_coarse_chan_bytes: gpubox_info.hdu_size * 4,
             num_timestep_coarse_chan_floats: gpubox_info.hdu_size,
+            num_timestep_coarse_chan_weight_floats: weight_floats,
             num_gpubox_files: gpubox_filenames.len(),
             legacy_conversion_table,
         })
@@ -472,6 +478,40 @@ impl CorrelatorContext {
         Ok(return_buffer)
     }
 
+    /// Read weights for a single timestep for a single coarse channel
+    /// The output weights are in order:
+    /// baseline,pol
+    ///
+    /// # Arguments
+    ///
+    /// * `corr_timestep_index` - index within the CorrelatorContext timestep array for the desired timestep. This corresponds
+    ///                      to the element within mwalibContext.timesteps.
+    ///
+    /// * `corr_coarse_chan_index` - index within the CorrelatorContext coarse_chan array for the desired coarse channel. This corresponds
+    ///                      to the element within mwalibContext.coarse_chans.
+    ///
+    ///
+    /// # Returns
+    ///
+    /// * A Result containing vector of 32 bit floats containing the data in [baseline][pol] order, if Ok.
+    ///
+    ///
+    pub fn read_weights_by_baseline(
+        &self,
+        corr_timestep_index: usize,
+        corr_coarse_chan_index: usize,
+    ) -> Result<Vec<f32>, GpuboxError> {
+        let mut return_buffer: Vec<f32> = vec![0.; self.num_timestep_coarse_chan_weight_floats];
+
+        self.read_weights_by_baseline_into_buffer(
+            corr_timestep_index,
+            corr_coarse_chan_index,
+            &mut return_buffer,
+        )?;
+
+        Ok(return_buffer)
+    }
+
     /// Validate input timestep_index and coarse_chan_index and return the fits_filename, batch index and hdu of the corresponding data
     ///
     /// # Arguments
@@ -678,6 +718,52 @@ impl CorrelatorContext {
         }
     }
 
+    /// Read weights from a single timestep for a single coarse channel
+    /// The output weights are in order:
+    /// baseline,pol
+    ///
+    /// # Arguments
+    ///
+    /// * `corr_timestep_index` - index within the CorrelatorContext timestep array for the desired timestep.
+    ///
+    /// * `corr_coarse_chan_index` - index within the CorrelatorContext coarse_chan array for the desired coarse channel.
+    ///
+    /// * `buffer` - Float buffer as a slice which will be filled with data from the HDU read in [baseline][pol] order.
+    ///
+    /// # Returns
+    ///
+    /// * A Result of Ok if success or a GpuboxError on failure.
+    ///
+    pub fn read_weights_by_baseline_into_buffer(
+        &self,
+        corr_timestep_index: usize,
+        corr_coarse_chan_index: usize,
+        buffer: &mut [f32],
+    ) -> Result<(), GpuboxError> {
+        // Validate input timestep_index and coarse_chan_index and return the fits_filename, batch index and hdu of the corresponding data
+        let (fits_filename, _, hdu_index) =
+            self.get_fits_filename_and_batch_and_hdu(corr_timestep_index, corr_coarse_chan_index)?;
+
+        // If we are not MWAXv2, just return an array of 1's
+        if self.mwa_version == MWAVersion::CorrMWAXv2 {
+            // Open the fits file
+            let mut fptr = fits_open!(&fits_filename)?;
+
+            // Use hdu_index + 1 as weights are always +1 from the visibilities for the same timestep
+            let hdu = fits_open_hdu!(&mut fptr, hdu_index + 1)?;
+
+            // Read into caller's buffer
+            get_fits_float_image_into_buffer!(&mut fptr, &hdu, buffer)?;
+
+            Ok(())
+        } else {
+            // Return an array of 1's
+            buffer.fill(1.0);
+
+            Ok(())
+        }
+    }
+
     /// Validates the first HDU of a gpubox file against metafits metadata
     ///
     /// In this case we call `validate_hdu_axes()`

src/correlator_context/test.rs

@@ -267,6 +267,12 @@ fn test_context_mwax() {
 
     assert_eq!(context.metafits_context.metafits_coarse_chans.len(), 24);
     assert_eq!(context.metafits_context.metafits_timesteps.len(), 120);
+
+    // weights is baselines x 4 visibility pols
+    assert_eq!(
+        context.num_timestep_coarse_chan_weight_floats,
+        ((128 * 129) / 2) * 4
+    )
 }
 
 #[test]
@@ -1181,3 +1187,123 @@ fn test_context_legacy_v1_get_fine_chan_feqs_some_coarse_chans() {
         fine_chan_freqs[128]
     );
 }
+
+#[test]
+fn test_read_weights_by_baseline_invalid_inputs() {
+    let mwax_metafits_filename = "test_files/1244973688_1_timestep/1244973688.metafits";
+    let mwax_filename = "test_files/1244973688_1_timestep/1244973688_20190619100110_ch114_000.fits";
+
+    // Open a context and load in a test metafits and gpubox file
+    let gpuboxfiles = vec![mwax_filename];
+    let context = CorrelatorContext::new(mwax_metafits_filename, &gpuboxfiles)
+        .expect("Failed to create CorrelatorContext");
+
+    // 99999 is invalid as a timestep for this observation
+    let result_invalid_timestep = context.read_by_baseline(99999, 10);
+
+    assert!(matches!(
+        result_invalid_timestep.unwrap_err(),
+        GpuboxError::InvalidTimeStepIndex(_)
+    ));
+
+    // 99999 is invalid as a coarse_chan for this observation
+    let result_invalid_coarse_chan = context.read_by_baseline(0, 99999);
+
+    assert!(matches!(
+        result_invalid_coarse_chan.unwrap_err(),
+        GpuboxError::InvalidCoarseChanIndex(_)
+    ));
+
+    // Do another test, this time with no data for the timestep/cc combo
+    let result = context.read_weights_by_baseline(10, 0);
+    assert!(result.is_err());
+    let error = result.unwrap_err();
+    assert!(
+        matches!(
+            error,
+            GpuboxError::NoDataForTimeStepCoarseChannel {
+                timestep_index: _,
+                coarse_chan_index: _
+            }
+        ),
+        "Error was {:?}",
+        error
+    );
+}
+
+#[test]
+fn test_read_weights_by_baseline_into_buffer_legacy() {
+    // Buffer should end up with 1.0's since legacy does not have ewights HDUs. Only MWAX does.
+    let mwax_metafits_filename = "test_files/1101503312_1_timestep/1101503312.metafits";
+    let mwax_filename =
+        "test_files/1101503312_1_timestep/1101503312_20141201210818_gpubox01_00.fits";
+
+    //
+    // Read the legacy file by frequency using mwalib
+    //
+    // Open a context and load in a test metafits and gpubox file
+    let gpuboxfiles = vec![mwax_filename];
+    let context = CorrelatorContext::new(mwax_metafits_filename, &gpuboxfiles)
+        .expect("Failed to create CorrelatorContext");
+
+    // Read into buffer by baseline
+    let mut mwalib_hdu_weights_by_bl1: Vec<f32> =
+        vec![0.; context.num_timestep_coarse_chan_weight_floats];
+
+    let result_read_bl_buffer =
+        context.read_weights_by_baseline_into_buffer(0, 0, &mut mwalib_hdu_weights_by_bl1);
+    assert!(result_read_bl_buffer.is_ok());
+
+    // Read by baseline
+    let mwalib_hdu_weights_by_bl2_result = context.read_weights_by_baseline(0, 0);
+    assert!(mwalib_hdu_weights_by_bl2_result.is_ok());
+    // Unwrap to a vector
+    let mwalib_hdu_weights_by_bl2 = mwalib_hdu_weights_by_bl2_result.unwrap();
+
+    // Expected result
+    let expected_weights_by_bl: Vec<f32> =
+        vec![1.0; context.num_timestep_coarse_chan_weight_floats];
+
+    // Check they all match each other
+    assert_eq!(mwalib_hdu_weights_by_bl1, expected_weights_by_bl);
+    assert_eq!(mwalib_hdu_weights_by_bl2, expected_weights_by_bl);
+}
+
+#[test]
+fn test_read_weights_by_baseline_into_buffer_mwax() {
+    // When this obs was created weights would all be 1.0 (but they are real unlike previous test for legacy MWA!)
+    let mwax_metafits_filename = "test_files/1244973688_1_timestep/1244973688.metafits";
+    let mwax_filename = "test_files/1244973688_1_timestep/1244973688_20190619100110_ch114_000.fits";
+
+    //
+    // Read the legacy file by frequency using mwalib
+    //
+    // Open a context and load in a test metafits and gpubox file
+    let gpuboxfiles = vec![mwax_filename];
+    let context = CorrelatorContext::new(mwax_metafits_filename, &gpuboxfiles)
+        .expect("Failed to create CorrelatorContext");
+
+    assert_eq!(8256 * 4, context.num_timestep_coarse_chan_weight_floats);
+
+    // Read into buffer by baseline
+    let mut mwalib_hdu_weights_by_bl1: Vec<f32> =
+        vec![0.; context.num_timestep_coarse_chan_weight_floats];
+
+    let result_read_bl_buffer =
+        context.read_weights_by_baseline_into_buffer(0, 10, &mut mwalib_hdu_weights_by_bl1);
+    assert!(result_read_bl_buffer.is_ok());
+
+    // Read by baseline
+    let mwalib_hdu_weights_by_bl2_result = context.read_weights_by_baseline(0, 10);
+    assert!(mwalib_hdu_weights_by_bl2_result.is_ok());
+    // Unwrap to a vector
+    let mwalib_hdu_weights_by_bl2 = mwalib_hdu_weights_by_bl2_result.unwrap();
+
+    // Expected result
+    let expected_weights_by_bl: Vec<f32> =
+        vec![1.0; context.num_timestep_coarse_chan_weight_floats];
+
+    // Check they all match each other
+    assert_eq!(mwalib_hdu_weights_by_bl1, expected_weights_by_bl);
+    assert_eq!(mwalib_hdu_weights_by_bl2, expected_weights_by_bl);
+}