Convert Matrix to a strided vector

tdms/include/arrays.h

@@ -229,75 +229,6 @@ struct DispersiveMultiLayer {
   bool is_dispersive(double near_zero_tolerance = 1e-15) const;
 };
 
-template<typename T>
-class Matrix {
-protected:
-  int n_rows = 0;
-  int n_cols = 0;
-  T **matrix = nullptr;
-
-public:
-  /**
-   * @brief Construct a new Matrix object, without assigned elements
-   *
-   */
-  Matrix() = default;
-  /**
-   * @brief Construct a new Matrix object, providing the dimensions
-   *
-   * @param n_rows,n_cols Number of rows and columns in the matrix
-   * @param initial_value The initial value of the elements, defaults to 0 to
-   * avoid initalised but unassigned values
-   */
-  Matrix(int n_rows, int n_cols) { allocate(n_rows, n_cols); }
-
-  inline T *operator[](int value) const { return matrix[value]; }
-  /**
-   * @brief Check whether this matrix has elements assigned
-   *
-   * @return true If this matrix has assigned elements
-   * @return false This matrix is currently unassigned
-   */
-  bool has_elements() { return matrix != nullptr; };
-
-  /**
-   * Allocate the memory for this matrix
-   *
-   * @param n_rows Number of rows
-   * @param n_cols Number of columns
-   */
-  void allocate(int n_rows, int n_cols, T initial_value = 0) {
-    this->n_rows = n_rows;
-    this->n_cols = n_cols;
-
-    matrix = (T **) malloc(sizeof(T *) * n_rows);
-    for (int i = 0; i < n_rows; i++) {
-      matrix[i] = (T *) malloc(sizeof(T) * n_cols);
-    }
-  };
-
-  int get_n_cols() const { return n_cols; }
-  int get_n_rows() const { return n_rows; }
-
-  /**
-   * Destructor. Must be defined in the header
-   */
-  ~Matrix() {
-    if (has_elements()) {
-      for (int i = 0; i < n_rows; i++) { free(matrix[i]); }
-      free(matrix);
-    }
-  };
-};
-
-class GratingStructure : public Matrix<int> {
-
-public:
-  GratingStructure(const mxArray *ptr, int I_tot);
-
-  ~GratingStructure();
-};
-
 template<typename T>
 class Vector {
 protected:
@@ -334,27 +265,6 @@ struct FrequencyVectors {
   std::vector<double> y;
 };
 
-/**
- * @brief Defines the numerical aperture of the objective, assuming that the
- * lens is centred on the origin of the PSTD simulation.
- *
- * In particular, since the fibre modes are imaged onto a Fourier plane of both
- * the physical fibre and the sample, the field scattered by the sample and
- * collected by the objective lens can have only a finite spatial support in the
- * aperature of the objective lens.
- *
- * Pupil[j][i] thus takes the value 1 for those (i,j) indices (note the order
- * swapping) within the aperture of the lens.
- */
-class Pupil : public Matrix<double> {
-public:
-  Pupil() = default;
-
-  void initialise(const mxArray *ptr, int n_rows, int n_cols);
-
-  ~Pupil();
-};
-
 /**
  * List of field components as integers
  */
@@ -373,20 +283,6 @@ class FieldComponentsVector : public Vector<int> {
   int index(int value);
 };
 
-class Vertices : public Matrix<int> {
-public:
-  Vertices() = default;
-
-  void initialise(const mxArray *ptr);
-
-  int n_vertices() { return n_rows; }
-
-  ~Vertices() {
-    if (has_elements()) { free_cast_matlab_2D_array(matrix); }
-    matrix = nullptr;
-  };
-};
-
 class DetectorSensitivityArrays {
 public:
   fftw_complex *v = nullptr;          // Flat fftw vector
@@ -398,11 +294,6 @@ class DetectorSensitivityArrays {
   ~DetectorSensitivityArrays();
 };
 
-/**
- * Matrix of c coefficients. See the pdf documentation for their definition
- */
-class CCoefficientMatrix : public Matrix<double> {};
-
 /**
  * Container for storing snapshots of the full-field
  */

tdms/include/arrays/tdms_matrix.h

@@ -0,0 +1,125 @@
+#pragma once
+
+#include <vector>
+
+#include "cell_coordinate.h"
+#include "matlabio.h"
+#include "matrix.h"
+
+/**
+ * @brief Template class for storing two dimensional data as a strided vector.
+ * @details Two dimensional data is stored as a strided vector, in row-major (C)
+ * format. The column dimension has the fastest-varying index, so
+ *
+ * this->(i, j) = data_[i * n_cols_ + j].
+ *
+ * @seealso This is consistent with the way hdf5 files store array-like data,
+ * see https://support.hdfgroup.org/HDF5/doc1.6/UG/12_Dataspaces.html.
+ * @tparam T Numerical datatype
+ */
+template<typename T>
+class Matrix {
+  friend class HDF5Reader;
+  friend class HDF5Writer;
+
+protected:
+  int n_rows_ = 0;
+  int n_cols_ = 0;
+
+  /*! Strided vector that will store the array data */
+  std::vector<T> data_;
+
+  /*! The total number of elements, according to the dimension values currently
+   * set. */
+  int n_elements() const { return n_rows_ * n_cols_; }
+
+public:
+  Matrix() = default;
+  Matrix(int n_rows, int n_cols) { allocate(n_rows, n_cols); }
+
+  int get_n_rows() const { return n_rows_; }
+  int get_n_cols() const { return n_cols_; }
+
+  /** @brief Subscript operator for the Matrix, retrieving the (i,j)-th element
+   */
+  T &operator()(int row, int col) { return data_[row * n_cols_ + col]; }
+  /** @copydoc operator() */
+  T operator()(int row, int col) const { return data_[row * n_cols_ + col]; }
+
+  /**
+   * @brief Allocate memory for this Matrix given the dimensions passed.
+   *
+   * @param n_rows,n_cols Number of rows and columns to assign.
+   */
+  void allocate(int n_rows, int n_cols) {
+    n_rows_ = n_rows;
+    n_cols_ = n_cols;
+    data_.resize(n_elements());
+  }
+
+  /**
+   * @brief Initialise this Matrix from a 2D-buffer of matching size.
+   * @details Data values are copied so that membership over this->data_ is
+   * preserved.
+   * @param buffer 2D buffer to read from
+   * @param n_rows,n_cols "Shape" of the read buffer to assign to this Matrix
+   * @param buffer_leads_n_cols If true, the buffer to read from is assumed to
+   * have dimensions [n_cols][n_rows]. If false, it is assumed to have
+   * dimensions [n_rows][n_cols].
+   */
+  void initialise(T **buffer, int n_rows, int n_cols,
+                  bool buffer_leads_n_cols = false) {
+    allocate(n_rows, n_cols);
+    if (buffer_leads_n_cols) {
+      for (int i = 0; i < n_rows; i++) {
+        for (int j = 0; j < n_cols; j++) { operator()(i, j) = buffer[j][i]; }
+      }
+    } else {
+      for (int i = 0; i < n_rows; i++) {
+        for (int j = 0; j < n_cols; j++) { operator()(i, j) = buffer[i][j]; }
+      }
+    }
+  }
+
+  /** @brief Whether this Matrix contains any elements */
+  bool has_elements() const { return n_elements() != 0; }
+};
+
+/** @brief Matrix of C-coefficients. See the pdf documentation for their
+ * definition. */
+typedef Matrix<double> CCoefficientMatrix;
+
+/**
+ * @brief Defines the numerical aperture of the objective, assuming that the
+ *        lens is centred on the origin of the PSTD simulation.
+ * @details In particular, since the fibre modes are imaged onto a Fourier plane
+ *         of both the physical fibre and the sample, the field scattered by the
+ * sample and collected by the objective lens can have only a finite spatial
+ * support in the aperture of the objective lens.
+ *
+ * Pupil(i, j) thus takes the value 1 for those (i,j) indices within the
+ * aperture of the lens.
+ */
+typedef Matrix<double> Pupil;
+
+/** TODO: Docstring */
+struct GratingStructure : public Matrix<int> {
+  GratingStructure(const mxArray *ptr, int I_tot);
+};
+
+/**
+ * @brief n_vertices-by-3 storage used for holding Yee cell indices.
+ * @details Each row of a Vertices instance consists of three integers (i, j, k)
+ * that form the index of a particular Yee cell in the simulation space.
+ */
+struct Vertices : public Matrix<int> {
+  Vertices() = default;
+
+  void initialise_from_matlab(const mxArray *ptr);
+
+  int n_vertices() { return n_rows_; }
+
+  /** @brief Convert the (i,j,k) index stored across the columns of row to an
+   * ijk struct. */
+  ijk index_in_row(int row) const;
+};

tdms/include/hdf5_io/hdf5_reader.h

@@ -5,6 +5,7 @@
 
 #include "arrays.h"
 #include "arrays/cuboid.h"
+#include "arrays/tdms_matrix.h"
 #include "interface.h"
 
 /**
@@ -13,7 +14,6 @@
  *          **double, but can be anything in general).
  */
 class HDF5Reader : public HDF5Base {
-
 public:
   /**
    * @brief Construct a new HDF5Reader for a named file.
@@ -81,19 +81,9 @@ class HDF5Reader : public HDF5Base {
               "Cannot read " + dataset_name + " into a 2D matrix, it has " +
               std::to_string(dimensions.size()) + " dimensions");
     }
-    int n_rows = dimensions[0];
-    int n_cols = dimensions[1];
-
-    SPDLOG_DEBUG("n_rows = {}; n_cols = {}", n_rows, n_cols);
-    T *buff = (T *) malloc(n_rows * n_cols * sizeof(T));
-    read(dataset_name, buff);
 
-    data_location.allocate(n_rows, n_cols);
-    for (unsigned int i = 0; i < n_rows; i++) {
-      for (unsigned int j = 0; j < n_cols; j++) {
-        data_location[i][j] = buff[i * n_cols + j];
-      }
-    }
+    data_location.allocate(dimensions[0], dimensions[1]);
+    read(dataset_name, data_location.data_.data());
     return;
   }
 

tdms/include/hdf5_io/hdf5_writer.h

@@ -2,7 +2,7 @@
 
 #include "hdf5_io/hdf5_base.h"
 
-#include "arrays.h"
+#include "arrays/tdms_matrix.h"
 
 class HDF5Writer : public HDF5Base {
 
@@ -22,7 +22,7 @@ class HDF5Writer : public HDF5Base {
    * @param size The size of the data array.
    * @param dimensions The number of dimensions of the array.
    */
-  void write(const std::string &dataname, double *data, int size,
+  void write(const std::string &dataname, const double *data, int size,
              hsize_t *dimensions);
 
   /**
@@ -35,16 +35,9 @@ class HDF5Writer : public HDF5Base {
    */
   template<typename T>
   void write(const std::string &dataname, const Matrix<T> &data) {
-    int n_cols = data.get_n_cols();
-    int n_rows = data.get_n_rows();
-    hsize_t dimension[2] = {static_cast<hsize_t>(n_rows),
-                            static_cast<hsize_t>(n_cols)};
-    T *buff = (T *) malloc(n_rows * n_cols * sizeof(T));
-    for (unsigned int i = 0; i < n_rows; i++) {
-      for (unsigned int j = 0; j < n_cols; j++) {
-        buff[i * n_cols + j] = data[i][j];
-      }
-    }
-    write(dataname, buff, 2, dimension);
+    hsize_t dimensions[2];
+    dimensions[0] = data.get_n_rows();
+    dimensions[1] = data.get_n_cols();
+    write(dataname, data.data_.data(), 2, dimensions);
   }
 };

tdms/include/simulation_manager/objects_from_infile.h

@@ -13,6 +13,7 @@
 #include "arrays/cuboid.h"
 #include "arrays/dtilde.h"
 #include "arrays/incident_field.h"
+#include "arrays/tdms_matrix.h"
 #include "cell_coordinate.h"
 #include "field.h"
 #include "fieldsample.h"
@@ -64,12 +65,12 @@ class IndependentObjectsFromInfile {
 
   IncidentField Ei;//< time-domain field
 
-  FrequencyVectors f_vec;                //< frequency vector
-  Pupil pupil;                           //< TODO
-  DTilde D_tilde;                        //< TODO
-  TDFieldExporter2D ex_td_field_exporter;//< two-dimensional field exporter
+  FrequencyVectors f_vec;//!< Frequency vector
+  Pupil pupil;           //!< Numerical aperture of the objective lens
+  DTilde D_tilde;        //!< TODO
+  TDFieldExporter2D ex_td_field_exporter;//!< two-dimensional field exporter
 
-  GridLabels input_grid_labels;//< cartesian labels of the Yee cells
+  GridLabels input_grid_labels;//!< cartesian labels of the Yee cells
 
   /* DERIVED VARIABLES FROM INDEPENDENT INPUTS */
 
@@ -90,7 +91,7 @@ class IndependentObjectsFromInfile {
  * inputs from this file.
  *
  * The Sources, GratingStructure, and FrequencyExtractVector can only be
- * initalised after the other arrays in the input file have been parsed.
+ * initialised after the other arrays in the input file have been parsed.
  *
  * As such, this object inherits the setup of IndependentObjectsFromInfile, and
  * then constructs the aforementioned arrays using a combination of information

tdms/include/simulation_manager/pstd_variables.h

@@ -7,7 +7,7 @@
 
 #include <fftw3.h>
 
-#include "arrays.h"
+#include "arrays/tdms_matrix.h"
 #include "cell_coordinate.h"
 
 /**

tdms/include/vertex_phasors.h

@@ -7,7 +7,7 @@
 
 #include <complex>
 
-#include "arrays.h"
+#include "arrays/tdms_matrix.h"
 #include "field.h"
 #include "grid_labels.h"