DispersiveMultiLayer no longer uses MATLAB (#286)

* Rework DispersiveMultiLayer class into a struct. TDMS builds. - xyz_vector struct of 3 vector<double>s introduced - DispersiveMultiLayer is now a struct with xyz_vector members - Preserve method for checking if the medium is dispersive - Temporarily disable unit tests for DispersiveMultiLayer class - HDF5Reader::read() can now assemble a DispersiveMultiLayer - HDF5Reader has a method for reading directly into a vector<T> - H5Dimension can be passed a H5::DataSet as well as a H5::DataSpace * Update paths now we need HDF5 data as well as matlab data * TDMS builds and tests HDF5Reader::read_data_from_group (passes, no) * Fix seg-fault in read_dataset_from_group * Add DispersiveMultiLayer test * Update DispersiveMultiLayer tests to not be entangled with the MATLAB jargon * Update the CI so things actually work * Move content of benchmark_scripts readme into doc/developers * Update CI to produce unit test `.mat` and `.hdf5` data (#289) * Move content of benchmark_scripts readme into doc/developers * Update CI attempt 1 * Actually use the right syntax * Force python version to prevent MacOS using python2.7 * Rename files to better match the data they produce * Aplease windows syntax * Remove pip cache due to known windows bug: actions/setup-python#436 * Apply suggestions from code review
UCL · May 19, 2023 · 30899d5 · 30899d5
1 parent eac765d
commit 30899d5
Show file tree

Hide file tree

Showing 27 changed files with 367 additions and 196 deletions.
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -56,10 +56,10 @@ jobs:
       - name: Set up MATLAB
         uses: matlab-actions/[email protected]
 
-      - name: Produce MATLAB unit test data
-        uses: matlab-actions/run-command@v1
+      - name: Setup Python
+        uses: actions/setup-python@v4
         with:
-          command: cd('tdms/tests/unit/matlab_benchmark_scripts/'), setup_unit_tests
+          python-version: '3.9'
 
       # -------------------------------------------------------------------------------
       # Ubuntu
@@ -128,6 +128,20 @@ jobs:
         run: |
           cmake --build . --config ${{ matrix.build_type }}
 
+      # -------------------------------------------------------------------------------
+      # Unit tests
+      - name: Produce MATLAB unit test data
+        if: matrix.build_testing == 'ON'
+        uses: matlab-actions/run-command@v1
+        with:
+          command: cd('tdms/tests/unit/benchmark_scripts/'), setup_unit_tests
+
+      - name: Produce hdf5 unit test data
+        if: matrix.build_testing == 'ON'
+        run: |
+          pip install -r ${{ github.workspace }}/tdms/tests/requirements.txt
+          python ${{ github.workspace }}/tdms/tests/unit/benchmark_scripts/create_hdf5_test_file.py
+
       - name: Run TDMS unit tests
         if: matrix.build_testing == 'ON'
         working-directory: ${{ runner.workspace }}/build
@@ -153,6 +167,8 @@ jobs:
         if: matrix.build_testing == 'ON'
         uses: codecov/codecov-action@v3
 
+      # -------------------------------------------------------------------------------
+      # Upload build artefact for system tests
       - name: Tar the build result to maintain permissions
         # https://github.com/actions/upload-artifact#maintaining-file-permissions-and-case-sensitive-files
         # https://github.com/actions/upload-artifact/issues/38

diff --git a/.gitignore b/.gitignore
@@ -15,6 +15,7 @@ html/
 **/.pytest_cache/
 tdms/tests/system/data/*.zip
 **.mat
+tdms/tests/unit/benchmark_scripts/unit_test_data/**
 
 # text editor files
 .idea/

diff --git a/doc/developers.md b/doc/developers.md
@@ -254,7 +254,22 @@ The doxygen-style comments will be included in this developer documentation.
 
 To run the unit tests, [compile](#compiling) with `-DBUILD_TESTING=ON`. Then run `ctest` from the build directory or execute the test executable `./tdms_tests`.
 
-It's good practice, and reassuring for your pull-request reviewers, if new C++ functionality is at covered by unit tests.
+It's good practice, and reassuring for your pull-request reviewers, if new C++ functionality is covered by unit tests.
+
+#### Benchmark Scripts and Data Generation
+
+The `tdms/tests/unit/benchmarking` directory contains scripts that produce input data for the unit tests, or that provide benchmarks for the units that are tested.
+
+The `C++` unit tests require the presence of a `.mat` or `.hdf5` file to read/write data from/to during testing.
+The locations of these files are coded into `tests/include/unit_test_utils.h`, but the files themselves are not committed to the repository - they can be created by running the `setup_unit_tests.m` and `create_hdf5_data.py` scripts.
+These scripts can then be updated to add/remove/edit the data available to the unit tests:
+- `create_hdf5_test_file.py` : Produces `hdf5_test_file.hdf5`; used when testing read/writing from `.hdf5` files.
+- `create_structure_array.m` : Produces `structure_array.mat`; used when testing reading/writing MATLAB `struct` arrays.
+- `create_class_data.m` : Produces `class_data.mat`; used when testing reading/writing `tdms` objects to/from `.mat` files.
+- `create_bad_class_data.m` : Produces `bad_class_data.mat`; used for testing failure cases when reading/writing `tdms` objects to/from `.mat` files.
+
+The `benchmark_` scripts perform band-limited interpolation (BLi) using `MATLAB`'s `interp` function.
+`TDMS`'s interpolation schemes are based off this `MATLAB` function (specficially, in the coefficients the scheme uses to interpolate), and are thus used to benchmark the accuracy of the scheme.
 
 ### Test coverage {#coverage}
 

diff --git a/tdms/include/arrays.h b/tdms/include/arrays.h
@@ -15,6 +15,13 @@
 #include "matlabio.h"
 #include "utils.h"
 
+template<typename T>
+struct xyz_vector {
+  std::vector<T> x = {};
+  std::vector<T> y = {};
+  std::vector<T> z = {};
+};
+
 template<typename T>
 class XYZTensor3D {
 public:
@@ -202,28 +209,23 @@ class DMaterial : public DCollectionBase, MaterialCollection {
   explicit DMaterial(const mxArray *ptr);
 };
 
-class DispersiveMultiLayer {
-public:
-  double *alpha = nullptr;
-  double *beta = nullptr;
-  double *gamma = nullptr;
-  XYZVectors kappa;
-  XYZVectors sigma;
-
+struct DispersiveMultiLayer {
 public:
-  explicit DispersiveMultiLayer(const mxArray *ptr);
+  std::vector<double> alpha;
+  std::vector<double> beta;
+  std::vector<double> gamma;
+  xyz_vector<double> kappa;
+  xyz_vector<double> sigma;
 
   /**
    * @brief Determines whether the (background) medium is dispersive
    *
-   * @param K_tot Number of Yee cells in the z-direction (number of entries in
-   * this->gamma)
    * @param near_zero_tolerance Tolerance for non-zero gamma (attenuation)
    * values
    * @return true Background is dispersive
    * @return false Background is not dispersive
    */
-  bool is_dispersive(int K_tot, double near_zero_tolerance = 1e-15);
+  bool is_dispersive(double near_zero_tolerance = 1e-15) const;
 };
 
 template<typename T>

diff --git a/tdms/include/hdf5_io/hdf5_dimension.h b/tdms/include/hdf5_io/hdf5_dimension.h
@@ -9,6 +9,7 @@ class H5Dimension : public std::vector<hsize_t> {
 public:
   H5Dimension() = default;
   H5Dimension(const H5::DataSpace &data_space);
+  H5Dimension(const H5::DataSet &data_set) : H5Dimension(data_set.getSpace()){};
 
   /**
    * @brief Whether these dimensions describe an array that is castable to a 1D
@@ -31,4 +32,11 @@ class H5Dimension : public std::vector<hsize_t> {
    * @return hsize_t
    */
   hsize_t max_dim() const { return *std::max_element(begin(), end()); };
+
+  /** @brief The total number of elements (product of the dimensions) */
+  int number_of_elements() const {
+    int product = 1;
+    for (hsize_t axis_size : *this) { product *= axis_size; }
+    return product;
+  }
 };
diff --git a/tdms/include/hdf5_io/hdf5_reader.h b/tdms/include/hdf5_io/hdf5_reader.h
@@ -23,9 +23,9 @@ class HDF5Reader : public HDF5Base {
       : HDF5Base(filename, H5F_ACC_RDONLY) {}
 
   /**
-   * @brief Read the dataset stored within a group into the buffer provided. Can
-   * be used to read MATLAB structs by treating the struct as the Group and
-   * field as the Dataset.
+   * @brief Read the dataset stored within a group into the buffer provided.
+   * @details Can be used to read MATLAB structs by treating the struct as the
+   * Group and field as the Dataset.
    * @tparam T C++ datatype to read data into.
    * @param group The Group within the file in which the dataset lives.
    * @param dataset The name of the dataset to fetch data from.
@@ -44,6 +44,34 @@ class HDF5Reader : public HDF5Base {
     requested_field.read(data, requested_field.getDataType());
   }
 
+  /**
+   * @brief Read the dataset stored within a group into the buffer provided,
+   * resizing the vector buffer accordingly.
+   * @details Can be used to read MATLAB structs by treating the struct as the
+   * Group and field as the Dataset.
+   * @tparam T C++ datatype to read data into.
+   * @param group The Group within the file in which the dataset lives.
+   * @param dataset The name of the dataset to fetch data from.
+   * @param[out] data The buffer into which to write the data.
+   */
+  template<typename T>
+  void read_dataset_in_group(const std::string &group,
+                             const std::string &dataset,
+                             std::vector<T> &data) const {
+    spdlog::debug("Reading {} from file: {}", group, filename_);
+
+    // Structs are saved as groups, so we need to fetch the group this struct is
+    // contained in
+    H5::Group structure_array = file_->openGroup(group);
+    // Then fetch the requested data and read it into the buffer provided,
+    // resizing the buffer if necessary
+    H5::DataSet requested_field = structure_array.openDataSet(dataset);
+    H5Dimension field_size(requested_field);
+    int number_of_elements = field_size.number_of_elements();
+    data.resize(number_of_elements);
+    requested_field.read(data.data(), requested_field.getDataType());
+  }
+
   /**
    * @brief Reads a named dataset from the HDF5 file.
    * @param dataname The name of the datset to be read.
@@ -139,4 +167,15 @@ class HDF5Reader : public HDF5Base {
    * @param cube
    */
   void read(Cuboid *cube) const;
+
+  /**
+   * @brief Read data from the file into a DispersiveMultiLayerObject
+   * @details Data is read from the "dispersive_aux" group. If this group does
+   * not exist, no data is written but no exception is thrown.
+   *
+   * If the group does exist, the alpha, beta, gamma, kappa, and sigma members
+   * are populated with the corresponding data entries.
+   * @param dml DispersiveMultiLayer object into which to write data.
+   */
+  void read(DispersiveMultiLayer *dml) const;
 };
diff --git a/tdms/src/arrays.cpp b/tdms/src/arrays.cpp
@@ -162,26 +162,9 @@ void DCollection::init_xyz_vectors(const mxArray *ptr, XYZVectors &arrays,
   }
 }
 
-DispersiveMultiLayer::DispersiveMultiLayer(const mxArray *ptr) {
-
-  if (mxIsEmpty(ptr)) { return; }
-  assert_is_struct_with_n_fields(ptr, 9, "dispersive_aux");
-
-  alpha = mxGetPr(ptr_to_vector_in(ptr, "alpha", "dispersive_aux"));
-  beta = mxGetPr(ptr_to_vector_in(ptr, "beta", "dispersive_aux"));
-  gamma = mxGetPr(ptr_to_vector_in(ptr, "gamma", "dispersive_aux"));
-  kappa.x = mxGetPr(ptr_to_matrix_in(ptr, "kappa_x", "dispersive_aux"));
-  kappa.y = mxGetPr(ptr_to_matrix_in(ptr, "kappa_y", "dispersive_aux"));
-  kappa.z = mxGetPr(ptr_to_matrix_in(ptr, "kappa_z", "dispersive_aux"));
-  sigma.x = mxGetPr(ptr_to_matrix_in(ptr, "sigma_x", "dispersive_aux"));
-  sigma.y = mxGetPr(ptr_to_matrix_in(ptr, "sigma_y", "dispersive_aux"));
-  sigma.z = mxGetPr(ptr_to_matrix_in(ptr, "sigma_z", "dispersive_aux"));
-}
-
-bool DispersiveMultiLayer::is_dispersive(int K_tot,
-                                         double near_zero_tolerance) {
-  for (int i = 0; i < K_tot; i++) {
-    if (fabs(gamma[i]) > near_zero_tolerance) {
+bool DispersiveMultiLayer::is_dispersive(double near_zero_tolerance) const {
+  for (double gamma_val : gamma) {
+    if (fabs(gamma_val) > near_zero_tolerance) {
       // non-zero attenuation constant of a Yee cell implies media is dispersive
       return true;
     }

diff --git a/tdms/src/hdf5_io/hdf5_reader.cpp b/tdms/src/hdf5_io/hdf5_reader.cpp
@@ -4,6 +4,8 @@
 #include <algorithm>
 #include <stdexcept>
 
+#include <spdlog/spdlog.h>
+
 using namespace std;
 
 void HDF5Reader::read(const string &plane, InterfaceComponent *ic) const {
@@ -52,3 +54,28 @@ void HDF5Reader::read(Cuboid *cube) const {
     cube->array[i] = (int) intermediate_buffer[i] - 1;
   }
 }
+
+void HDF5Reader::read(DispersiveMultiLayer *dml) const {
+  string group_name = "dispersive_aux";
+  // Deal with the case of an empty input
+  if (!file_->nameExists(group_name)) {
+    spdlog::info(group_name + " is not a group: assuming empty input");
+    return;
+  } else {
+    // This is a group - it should have 9 members and we can quickly check this
+    H5::Group dispersive_aux = file_->openGroup(group_name);
+    if (dispersive_aux.getNumObjs() != 9) {
+      throw runtime_error("dispersive_aux does not have exactly 9 members!");
+    }
+  }
+  // Assuming non-empty input, setup the data appropriately
+  read_dataset_in_group(group_name, "alpha", dml->alpha);
+  read_dataset_in_group(group_name, "beta", dml->beta);
+  read_dataset_in_group(group_name, "gamma", dml->gamma);
+  read_dataset_in_group(group_name, "kappa_x", dml->kappa.x);
+  read_dataset_in_group(group_name, "kappa_y", dml->kappa.y);
+  read_dataset_in_group(group_name, "kappa_z", dml->kappa.z);
+  read_dataset_in_group(group_name, "sigma_x", dml->sigma.x);
+  read_dataset_in_group(group_name, "sigma_y", dml->sigma.y);
+  read_dataset_in_group(group_name, "sigma_z", dml->sigma.z);
+}
diff --git a/tdms/src/simulation_manager/objects_from_infile.cpp b/tdms/src/simulation_manager/objects_from_infile.cpp
@@ -17,9 +17,7 @@ IndependentObjectsFromInfile::IndependentObjectsFromInfile(
       Dmaterial(matrices_from_input_file["Dmaterial"]),// get Dmaterial
       C(matrices_from_input_file["C"]),                // get C
       D(matrices_from_input_file["D"]),                // get D
-      matched_layer(
-              matrices_from_input_file["dispersive_aux"]),// get dispersive_aux
-      Ei(matrices_from_input_file["tdfield"])             // get tdfield
+      Ei(matrices_from_input_file["tdfield"])          // get tdfield
 {
   /* Set FDTD/PSTD-dependent variable skip_tdf [1: PSTD, 6: FDTD] */
   skip_tdf = in_flags["use_pstd"] ? 1 : 6;
@@ -140,9 +138,7 @@ IndependentObjectsFromInfile::IndependentObjectsFromInfile(
   }
 
   // work out if we have a dispersive background
-  if (params.is_disp_ml) {
-    params.is_disp_ml = matched_layer.is_dispersive(IJK_tot.k);
-  }
+  if (params.is_disp_ml) { params.is_disp_ml = matched_layer.is_dispersive(); }
 
   // Set dt so that an integer number of time periods fits within a sinusoidal
   // period

diff --git a/tdms/tests/include/array_test_class.h b/tdms/tests/include/array_test_class.h
@@ -60,24 +60,6 @@ class DetectorSensitivityArraysTest : public AbstractArrayTest {
   std::string get_class_name() override { return "DetectorSensitivityArray"; }
 };
 
-/** @brief Unit tests for DispersiveMultilayer */
-class DispersiveMultilayerTest : public AbstractArrayTest {
-private:
-  const int n_fields = 9;
-  const char *fieldnames[9] = {"alpha",   "beta",    "gamma",
-                               "kappa_x", "kappa_y", "kappa_z",
-                               "sigma_x", "sigma_y", "sigma_z"};
-
-  void test_empty_construction() override;
-  void test_wrong_input_type() override;
-  void test_correct_construction() override;
-  // test: is_dispersive()
-  void test_other_methods() override;
-
-public:
-  std::string get_class_name() override { return "DispersiveMultilayer"; }
-};
-
 // Test methods check the performance of initialise, as this is the de-facto
 // constructor
 class DTildeTest : public AbstractArrayTest {

diff --git a/tdms/tests/include/unit_test_utils.h b/tdms/tests/include/unit_test_utils.h
@@ -17,41 +17,30 @@ namespace tdms_unit_test_data {
 
 #ifdef CMAKE_SOURCE_DIR
 inline std::string tdms_object_data(std::string(CMAKE_SOURCE_DIR) +
-                                    "/tests/unit/matlab_benchmark_scripts/"
-                                    "matlab_data/class_data.mat");
+                                    "/tests/unit/benchmark_scripts/"
+                                    "unit_test_data/class_data.mat");
 inline std::string tdms_bad_object_data(std::string(CMAKE_SOURCE_DIR) +
-                                        "/tests/unit/matlab_benchmark_scripts/"
-                                        "matlab_data/bad_class_data.mat");
-#else
-inline std::string tdms_object_data(std::filesystem::current_path() /
-                                    "../tests/unit/matlab_benchmark_scripts/"
-                                    "matlab_data/class_data.mat");
-inline std::string
-        tdms_bad_object_data(std::filesystem::current_path() /
-                             "../tests/unit/matlab_benchmark_scripts/"
-                             "matlab_data/bad_class_data.mat");
-#endif
-
-/* The struct_testdata is a .mat file containing a single structure array,
-example_struct. This structure array has the following data.
-
-MEMBER: value, MATLAB type
-double_no_decimal: 1.0, double
-double_half: 0.5, double
-string: 'tdms', char array
-boolean: 1, logical
-uint_345: 3 * 4 * 5 array of 1s, uint8
-double_22: [0.25, 0.5; 0.75, 1.], double
-complex_22: [0., -i; i, 0], complex
-*/
-#ifdef CMAKE_SOURCE_DIR
+                                        "/tests/unit/benchmark_scripts/"
+                                        "unit_test_data/bad_class_data.mat");
+inline std::string hdf5_test_file(std::string(CMAKE_SOURCE_DIR) +
+                                  "/tests/unit/benchmark_scripts/"
+                                  "unit_test_data/hdf5_test_file.hdf5");
 inline std::string struct_testdata(std::string(CMAKE_SOURCE_DIR) +
-                                   "/tests/unit/matlab_benchmark_scripts/"
-                                   "matlab_data/structure_array.mat");
+                                   "/tests/unit/benchmark_scripts/"
+                                   "unit_test_data/structure_array.mat");
 #else
+inline std::string tdms_object_data(std::filesystem::current_path() /
+                                    "../tests/unit/benchmark_scripts/"
+                                    "unit_test_data/class_data.mat");
+inline std::string tdms_bad_object_data(std::filesystem::current_path() /
+                                        "../tests/unit/benchmark_scripts/"
+                                        "unit_test_data/bad_class_data.mat");
+inline std::string hdf5_test_file(std::filesystem::current_path() /
+                                  "../tests/unit/benchmark_scripts/"
+                                  "unit_test_data/hdf5_test_file.hdf5");
 inline std::string struct_testdata(std::filesystem::current_path() /
-                                   "../tests/unit/matlab_benchmark_scripts/"
-                                   "matlab_data/structure_array.mat");
+                                   "../tests/unit/benchmark_scripts/"
+                                   "unit_test_data/structure_array.mat");
 #endif
 
 }// namespace tdms_unit_test_data