API changes and small improvements.

The FFI interface has the biggest API changes. It is now allowing errors
to be propagated to the caller, and is overall a better C API. Examples
have been updated to demonstrate this new API.

The biggest change to the Rust API is returning a Vec instead of Array1
from the calc_jones*_array functions. I'm not sure why I wanted Array1
in the first place; Vec is much more useful here.

The Python API now allows the same calc_jones or calc_jones_array
functions to be used with either 16 or 32 dipole gains.

Move FEE tests into their own file. Also move the FFI code under the FEE
directory, as it's all FEE-related anyway.

Add a binary "verify-beam-file". It attempts to detect defects in a
user's HDF5 files by calculating beam responses, which trigger
calculating m_accums and n_accums. If there are problems in generating
those, errors should appear.

Add a function to "fix" dipole gains ("amps"). Tests are also included.

In calc_sigmas and calc_jones_direct, using the for_each iterator
method *could* be faster than using a for loop, so why not.

Using get_unchecked appears to remove boundary checks when accessing
J_POWER_TABLE; it's totally safe to do this and apparently makes things
a little faster. Elsewhere, using an iterator instead of accessing Jones
matrix elements directly might have the same effect.

Make the freqs member of FEEBeam private. This prevents unintentional
modifications and problems. Access to the freqs is given instead by a
get_freqs method.

Using to_bits on a float is better than converting a float to an int for
hashing purposes.

Small changes elsewhere, like tidying documentation.

Cargo.toml

@@ -34,7 +34,7 @@ hdf5 = "0.7.*"
 marlu = { version = "0.3.*", default-features = false }
 thiserror = "1.0.*"
 
-hdf5-sys = { version = "0.7.1", features = ["static"], optional = true }
+hdf5-sys = { version = "0.7.*", features = ["static", "threadsafe"], optional = true }
 
 pyo3 = { version = "0.13.*", features = ["extension-module"], optional = true }
 numpy = { version = "0.13.0", optional = true }

build.rs

@@ -10,7 +10,7 @@ fn main() {
  // Generate a C header for hyperbeam and write it to the include
  // directory. This routine only need to be done if the ffi module has
  // changed.
- println!("cargo:rerun-if-changed=src/ffi.rs");
+ println!("cargo:rerun-if-changed=src/fee/ffi/mod.rs");
  // Only do this if we're not on docs.rs (doesn't like writing files outside
  // of OUT_DIR).
  match env::var("DOCS_RS").as_deref() {

examples/beam_calcs.c

@@ -6,6 +6,7 @@
 // gcc -O3 -I ../include/ -L ../target/release/ -l mwa_hyperbeam ./beam_calcs.c -o beam_calcs
 // LD_LIBRARY_PATH=../target/release ./beam_calcs ../mwa_full_embedded_element_pattern.h5
 
+#include <complex.h>
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
@@ -18,8 +19,13 @@ int main(int argc, char *argv[]) {
  exit(1);
  }
 
- // Get a new beam from hyperbeam.
- FEEBeam *beam = new_fee_beam(argv[1]);
+ // Get a new beam object from hyperbeam.
+ FEEBeam *beam;
+ char error[200];
+ if (new_fee_beam(argv[1], &beam, error)) {
+ printf("Got an error when trying to make an FEEBeam: %s\n", error);
+ return EXIT_FAILURE;
+ }
 
  // Set up the direction and pointing to test.
  double az = 45.0 * M_PI / 180.0;
@@ -37,31 +43,40 @@ int main(int argc, char *argv[]) {
  // https://github.com/JLBLine/polarisation_tests_for_FEE
  int parallactic = 1;
 
- // Calculate the Jones matrix for this direction and pointing.
- double *jones = calc_jones(beam, az, za, freq_hz, delays, amps, norm_to_zenith, parallactic);
+ // Calculate the Jones matrix for this direction and pointing. This Jones
+ // matrix is on the stack.
+ complex double jones[4];
+ // hyperbeam expects a pointer to doubles. Casting the pointer works fine.
+ if (calc_jones(beam, az, za, freq_hz, delays, amps, 16, norm_to_zenith, parallactic, (double *)&jones, error)) {
+ printf("Got an error when running calc_jones: %s\n", error);
+ return EXIT_FAILURE;
+ }
  printf("The returned Jones matrix:\n");
- printf("[[%+.8f%+.8fi,", jones[0], jones[1]);
- printf(" %+.8f%+.8fi]\n", jones[2], jones[3]);
- printf(" [%+.8f%+.8fi,", jones[4], jones[5]);
- printf(" %+.8f%+.8fi]]\n", jones[6], jones[7]);
+ printf("[[%+.8f%+.8fi,", creal(jones[0]), cimag(jones[0]));
+ printf(" %+.8f%+.8fi]\n", creal(jones[1]), cimag(jones[1]));
+ printf(" [%+.8f%+.8fi,", creal(jones[2]), cimag(jones[2]));
+ printf(" %+.8f%+.8fi]]\n", creal(jones[3]), cimag(jones[3]));
 
- // Amps can have 32 elements to specify X and Y elements of the dipoles, but
- // another function is needed to use this extra info. The first 16 elements
- // are X elements, second 16 are Y elements.
+ // Amps can have 32 elements to specify X and Y elements of the dipoles. The
+ // first 16 elements are X elements, second 16 are Y elements.
  double amps_2[32] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1,
  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1};
- double *jones_2 = calc_jones_all_amps(beam, az, za, freq_hz, delays, amps_2, norm_to_zenith, parallactic);
+ // This Jones matrix is on the heap.
+ complex double *jones_2 = malloc(4 * sizeof(complex double));
+ if (calc_jones(beam, az, za, freq_hz, delays, amps_2, 32, norm_to_zenith, parallactic, (double *)jones_2, error)) {
+ printf("Got an error when running calc_jones_all_amps: %s\n", error);
+ return EXIT_FAILURE;
+ }
  // The resulting Jones matrix has different elements on the second row,
  // corresponding to the Y element; this is because we only altered the Y
  // amps.
  printf("The returned Jones matrix with altered Y amps:\n");
- printf("[[%+.8f%+.8fi,", jones_2[0], jones_2[1]);
- printf(" %+.8f%+.8fi]\n", jones_2[2], jones_2[3]);
- printf(" [%+.8f%+.8fi,", jones_2[4], jones_2[5]);
- printf(" %+.8f%+.8fi]]\n", jones_2[6], jones_2[7]);
+ printf("[[%+.8f%+.8fi,", creal(jones_2[0]), cimag(jones_2[0]));
+ printf(" %+.8f%+.8fi]\n", creal(jones_2[1]), cimag(jones_2[1]));
+ printf(" [%+.8f%+.8fi,", creal(jones_2[2]), cimag(jones_2[2]));
+ printf(" %+.8f%+.8fi]]\n", creal(jones_2[3]), cimag(jones_2[3]));
 
- // Free the Jones matrices.
- free(jones);
+ // Free the heap-allocated Jones matrix.
  free(jones_2);
  // Free the beam - we must use a special function to do this.
  free_fee_beam(beam);

examples/beam_calcs.py

@@ -46,11 +46,10 @@
 print(jones[0])
 
 # It's also possible to supply amps for all dipole elements. The first 16 amps
-# are for X elements, the second 16 are for Y elements. beam.calc_jones_all_amps
-# is also available.
-amps = [1.0] * 32
+# are for X elements, the second 16 are for Y elements.
+amps = np.ones(32)
 amps[-1] = 0
-jones = beam.calc_jones_array_all_amps(
+jones = beam.calc_jones_array(
  az[:1], za[:1], freq, delays, amps, norm_to_zenith, parallactic)
 print("First Jones matrix with altered Y amps:")
 print(jones[0])

examples/beam_calcs.rs

@@ -83,17 +83,15 @@ fn main() -> Result<(), anyhow::Error> {
  results
  };
  println!("The first Jones matrix:");
- for j in jones.iter() {
- // This works, but the formatting for this isn't very pretty.
- // println!("{}", j);
+ // This works, but the formatting for this isn't very pretty.
+ // println!("{}", jones[0]);
 
- // For demonstrations' sake, this gives easier-to-read output.
- println!(
- "[[{:+.8}, {:+.8}]\n [{:+.8}, {:+.8}]]",
- j[0], j[1], j[2], j[3]
- );
- break;
- }
+ // For demonstrations' sake, this gives easier-to-read output.
+ let j = jones[0];
+ println!(
+ "[[{:+.8}, {:+.8}]\n [{:+.8}, {:+.8}]]",
+ j[0], j[1], j[2], j[3]
+ );
 
  Ok(())
 }

examples/beam_calcs_many.c

@@ -8,6 +8,7 @@
 // gcc -O3 -I ../include/ -L ../target/release/ -l mwa_hyperbeam ./beam_calcs_many.c -o beam_calcs_many
 // LD_LIBRARY_PATH=../target/release ./beam_calcs_many ../mwa_full_embedded_element_pattern.h5
 
+#include <complex.h>
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
@@ -20,8 +21,13 @@ int main(int argc, char *argv[]) {
  exit(1);
  }
 
- // Get a new beam from hyperbeam.
- FEEBeam *beam = new_fee_beam(argv[1]);
+ // Get a new beam object from hyperbeam.
+ FEEBeam *beam;
+ char error[200];
+ if (new_fee_beam(argv[1], &beam, error)) {
+ printf("Got an error when trying to make an FEEBeam: %s\n", error);
+ return EXIT_FAILURE;
+ }
 
  // Set up the directions to test.
  int num_directions = 5000;
@@ -45,22 +51,32 @@ int main(int argc, char *argv[]) {
 
  // Calculate the Jones matrices for all directions. Rust will do this in
  // parallel.
- double *jones = calc_jones_array(beam, num_directions, az, za, freq_hz, delays, amps, norm_to_zenith, parallactic);
+ complex double *jones;
+ // hyperbeam expects a pointer to doubles. Casting the pointer works fine.
+ if (calc_jones_array(beam, num_directions, az, za, freq_hz, delays, amps, 16, norm_to_zenith, parallactic,
+ (double **)&jones, error)) {
+ printf("Got an error when running calc_jones_array: %s\n", error);
+ return EXIT_FAILURE;
+ }
  printf("The first Jones matrix:\n");
- printf("[[%+.8f%+.8fi,", jones[0], jones[1]);
- printf(" %+.8f%+.8fi]\n", jones[2], jones[3]);
- printf(" [%+.8f%+.8fi,", jones[4], jones[5]);
- printf(" %+.8f%+.8fi]]\n", jones[6], jones[7]);
+ printf("[[%+.8f%+.8fi,", creal(jones[0]), cimag(jones[0]));
+ printf(" %+.8f%+.8fi]\n", creal(jones[1]), cimag(jones[1]));
+ printf(" [%+.8f%+.8fi,", creal(jones[2]), cimag(jones[2]));
+ printf(" %+.8f%+.8fi]]\n", creal(jones[3]), cimag(jones[3]));
 
  double amps_2[32] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0};
- double *jones_2 =
- calc_jones_array_all_amps(beam, num_directions, az, za, freq_hz, delays, amps_2, norm_to_zenith, parallactic);
+ complex double *jones_2;
+ if (calc_jones_array(beam, num_directions, az, za, freq_hz, delays, amps_2, 32, norm_to_zenith, parallactic,
+ (double **)&jones_2, error)) {
+ printf("Got an error when running calc_jones_array_all_amps: %s\n", error);
+ return EXIT_FAILURE;
+ }
  printf("The first Jones matrix with altered Y amps:\n");
- printf("[[%+.8f%+.8fi,", jones_2[0], jones_2[1]);
- printf(" %+.8f%+.8fi]\n", jones_2[2], jones_2[3]);
- printf(" [%+.8f%+.8fi,", jones_2[4], jones_2[5]);
- printf(" %+.8f%+.8fi]]\n", jones_2[6], jones_2[7]);
+ printf("[[%+.8f%+.8fi,", creal(jones_2[0]), cimag(jones_2[0]));
+ printf(" %+.8f%+.8fi]\n", creal(jones_2[1]), cimag(jones_2[1]));
+ printf(" [%+.8f%+.8fi,", creal(jones_2[2]), cimag(jones_2[2]));
+ printf(" %+.8f%+.8fi]]\n", creal(jones_2[3]), cimag(jones_2[3]));
 
  // Freeing memory.
  free(az);

examples/beam_calcs_many_omp.c

@@ -9,6 +9,7 @@
 // gcc -O3 -fopenmp -I ../include/ -L ../target/release/ -l mwa_hyperbeam ./beam_calcs_many_omp.c -o beam_calcs_many_omp
 // LD_LIBRARY_PATH=../target/release ./beam_calcs_many_omp ../mwa_full_embedded_element_pattern.h5
 
+#include <complex.h>
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
@@ -21,8 +22,13 @@ int main(int argc, char *argv[]) {
  exit(1);
  }
 
- // Get a new beam from hyperbeam.
- FEEBeam *beam = new_fee_beam(argv[1]);
+ // Get a new beam object from hyperbeam.
+ FEEBeam *beam;
+ char error[200];
+ if (new_fee_beam(argv[1], &beam, error)) {
+ printf("Got an error when trying to make an FEEBeam: %s\n", error);
+ return EXIT_FAILURE;
+ }
 
  // Set up the directions to test.
  int num_directions = 5000;
@@ -47,24 +53,25 @@ int main(int argc, char *argv[]) {
  int parallactic = 1;
 
  // Calculate the Jones matrices for all directions.
- double **jones = malloc(num_directions * 8 * sizeof(double));
+ complex double *jones = malloc(num_directions * 4 * sizeof(complex double));
 #pragma omp parallel for
  for (int i = 0; i < num_directions; i++) {
- double *j = calc_jones(beam, az[i], za[i], freq_hz, delays, amps, norm_to_zenith, parallactic);
- jones[i] = j;
+ // hyperbeam expects a pointer to doubles. Casting the pointer works fine.
+ if (calc_jones(beam, az[i], za[i], freq_hz, delays, amps, 16, norm_to_zenith, parallactic,
+ (double *)(jones + i * 4), error)) {
+ printf("Got an error when running calc_jones: %s\n", error);
+ return EXIT_FAILURE;
+ }
  }
  printf("The first Jones matrix:\n");
- printf("[[%+.8f%+.8fi,", jones[0][0], jones[0][1]);
- printf(" %+.8f%+.8fi]\n", jones[0][2], jones[0][3]);
- printf(" [%+.8f%+.8fi,", jones[0][4], jones[0][5]);
- printf(" %+.8f%+.8fi]]\n", jones[0][6], jones[0][7]);
+ printf("[[%+.8f%+.8fi,", creal(jones[0]), cimag(jones[0]));
+ printf(" %+.8f%+.8fi]\n", creal(jones[1]), cimag(jones[1]));
+ printf(" [%+.8f%+.8fi,", creal(jones[2]), cimag(jones[2]));
+ printf(" %+.8f%+.8fi]]\n", creal(jones[3]), cimag(jones[3]));
 
  // Freeing memory.
  free(az);
  free(za);
- for (int i = 0; i < num_directions; i++) {
- free(jones[i]);
- }
  free(jones);
 
  // Free the beam - we must use a special function to do this.

examples/get_freqs.c

@@ -18,15 +18,20 @@ int main(int argc, char *argv[]) {
  exit(1);
  }
 
- // Get a new beam from hyperbeam.
- FEEBeam *beam = new_fee_beam(argv[1]);
+ // Get a new beam object from hyperbeam.
+ FEEBeam *beam;
+ char error[200];
+ if (new_fee_beam(argv[1], &beam, error)) {
+ printf("Got an error when trying to make an FEEBeam: %s\n", error);
+ }
 
  // Of the available frequencies, which is closest to 255 MHz?
  printf("Closest freq. to 255 MHz: %.2f MHz\n", (double)closest_freq(beam, 255000000) / 1e6);
 
  // Get the frequecies from the FEEBeam struct.
- unsigned num_freqs = get_num_fee_beam_freqs(beam);
- unsigned *freqs = get_fee_beam_freqs(beam);
+ size_t num_freqs;
+ const unsigned *freqs;
+ get_fee_beam_freqs(beam, &freqs, &num_freqs);
 
  // Print them out.
  printf("All frequencies:\n");
@@ -42,8 +47,7 @@ int main(int argc, char *argv[]) {
  }
  }
 
- // We own the freqs array - free it.
- free(freqs);
+ // We DON'T own the freqs array - don't free it.
  // Free the beam - we must use a special function to do this.
  free_fee_beam(beam);
 

src/bin/verify-beam-file.rs

@@ -0,0 +1,38 @@
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+//! This program calculates the "m_accum" and "n_accum" for all frequencies
+//! specified in the given HDF5 file to verify the file's data. hyperbeam used
+//! to have runtime checks for this, but it's much easier on the computer to not
+//! check something that likely never happens.
+
+use mwa_hyperbeam::fee::{FEEBeam, InitFEEBeamError};
+
+fn main() {
+ // Test each input file.
+ for beam_file in std::env::args().skip(1) {
+ // If this threw an error, it was during initialisation.
+ if let Err(e) = test_file(&beam_file) {
+ println!("File '{}' failed to create an FEEBeam: {}", &beam_file, e);
+ }
+ }
+}
+
+fn test_file(beam_file: &str) -> Result<(), InitFEEBeamError> {
+ println!("Testing file '{}'", beam_file);
+ let beam = FEEBeam::new(&beam_file)?;
+ // It does not matter what the direction, delays or amps are, the m_accum
+ // and n_accum that we're interested in testing depend only on the
+ // frequency. (Frequency determined the number of dipole coeffs, which
+ // determines which m_accum and n_accum values are used.)
+ for &file_freq in beam.get_freqs() {
+ println!("Testing freq {}", file_freq);
+ // If this blows up, we know there's a problem...
+ beam.calc_jones_eng(0.0, 0.0, file_freq, &[0; 16], &[1.0; 16], false)
+ .unwrap();
+ }
+
+ println!("File '{}' is all good!", beam_file);
+ Ok(())
+}

src/fee/error.rs

@@ -9,7 +9,10 @@ use thiserror::Error;
 #[derive(Error, Debug)]
 #[allow(clippy::upper_case_acronyms)]
 pub enum InitFEEBeamError {
- #[error("One of HDF5 datasets started with 'X_'; what's wrong with your file?")]
+ #[error("Specified beam file '{0}' doesn't exist")]
+ BeamFileDoesntExist(String),
+
+ #[error("One of the HDF5 datasets started with 'X_'; what's wrong with your file?")]
  MissingDipole,
 
  #[error("No HDF5 datasets started with a 'X'; is there any data in the file?")]
@@ -35,7 +38,7 @@ pub enum InitFEEBeamError {
 
  /// An error associated with the hdf5 crate.
  #[error("HDF5 error: {0}")]
- Hdf5Error(#[from] hdf5::Error),
+ Hdf5(#[from] hdf5::Error),
 }
 
 #[derive(Error, Debug)]