rev: split handlers; handle size invariant

CMakeLists.txt

@@ -231,10 +231,8 @@ add_subdirectory("src/realizations/catchment")
 add_subdirectory("src/models/kernels/reservoir")
 add_subdirectory("src/models/kernels/evapotranspiration")
 add_subdirectory("src/forcing")
-
-add_library(ngen_utils INTERFACE)
-add_library(NGen::utils ALIAS ngen_utils)
-target_include_directories(ngen_utils INTERFACE ${CMAKE_CURRENT_LIST_DIR}/include/utilities)
+add_subdirectory("src/utilities/mdarray")
+add_subdirectory("src/utilities/mdframe")
 
 target_link_libraries(ngen PUBLIC
         NGen::core

include/utilities/mdarray/mdarray.hpp

@@ -1,6 +1,7 @@
 #ifndef NGEN_MDARRAY_DEFINITION_HPP
 #define NGEN_MDARRAY_DEFINITION_HPP
 
+#include <type_traits>
 #include <vector>
 
 #include <span.hpp>
@@ -23,118 +24,125 @@ class mdarray
     struct iterator;
     friend iterator;
 
-    // Deleted default constructor. mdarray must be given a rank at least.
-    mdarray() = delete;
+  private:
+    template<typename Tp>
+    inline size_type max_size(Tp shape)
+    {
+        size_type max = 1;
+        for (size_type dimension : shape) {
+            max *= dimension;
+        }
+        return max;
+    }
+
+  public:
 
-    mdarray(size_type rank) noexcept
-        : m_shape(rank)
-        , m_data() {};
+    mdarray() = default;
 
-    mdarray(boost::span<size_t> dsizes)
+    mdarray(const boost::span<const size_type> dsizes)
         : m_shape(dsizes.begin(), dsizes.end())
-        , m_data(this->max_size()) {};
+        , m_data(this->max_size(dsizes)) {};
 
-    mdarray(std::initializer_list<size_t> dsizes)
+    constexpr mdarray(std::initializer_list<size_type> dsizes)
         : m_shape(dsizes.begin(), dsizes.end())
-        , m_data(this->max_size()) {};
+        , m_data(this->max_size(dsizes)) {};
 
     /**
      * Retrieve a reference to the value at the given index.
      *
      * @example
-     * mdarray<int> x(2); // rank 2
-     * x.at({0, 0});      // get the value at (0, 0) in a 2D array.
-     * x.at({1, 1}) = 3;  // set the value at (1, 1) in a 2D array.
+     * mdarray<int> x({2, 2}); // rank 2
+     * x.at({0, 0});           // get the value at (0, 0) in a 2D array.
+     * x.at({1, 1}) = 3;       // set the value at (1, 1) in a 2D array.
      *
      * @param n Index list
      * @return reference 
      */
-    reference at(const boost::span<const size_t> n)
+    reference at(const boost::span<const size_type> n)
     {
+        this->bounds_check(n);
+
         return this->m_data.at(this->index(n));
     }
 
     /**
-     * Retrieve const reference to the value at the given index.
+     * Retrieve a reference to the value at the given index.
      *
      * @example
-     * mdarray<int> x(2); // rank 2
-     * x.at({0, 0});      // get the value at (0, 0) in a 2D array.
+     * mdarray<int> x({2, 2}); // rank 2
+     * x.at({0, 0});           // get the value at (0, 0) in a 2D array.
+     * x.at({1, 1}) = 3;       // set the value at (1, 1) in a 2D array.
      *
      * @param n Index list
      * @return reference 
      */
-    const_reference at(const boost::span<const size_t> n) const
+    reference operator[](const boost::span<const size_type> n) noexcept
     {
-        return this->m_data.at(this->index(n));
+        return this->m_data[this->index(n)];
     }
 
-    void insert(const boost::span<const size_t> n, value_type value)
+    /**
+     * Retrieve const reference to the value at the given index.
+     *
+     * @example
+     * mdarray<int> x({1, 1}); // rank 2
+     * x.at({0, 0});           // get the value at (0, 0) in a 2D array.
+     *
+     * @param n Index list
+     * @return reference 
+     */
+    const_reference at(const boost::span<const size_type> n) const
     {
-        for (size_type i = 0; i < this->m_shape.size(); i++) {
-            if (n[i] >= this->m_shape[i])
-                throw std::out_of_range(
-                    "index " + std::to_string(n[i]) +
-                    " must be less than dimension size " +
-                    std::to_string(this->m_shape[i])
-                );
-        }
+        this->bounds_check(n);
 
-        size_type index = this->index(n);
-
-        if (index >= this->max_size())
-            throw std::out_of_range(
-                "index " + std::to_string(index) +
-                " is larger than maximum possible address index (" +
-                std::to_string(this->max_size()) + ")"
-            );
-
-        if (index >= this->size() || this->size() == 0)
-            this->m_data.resize(index + 1);
-
-        this->m_data[index] = value;
+        return this->m_data.at(this->index(n));
     }
 
-    void insert(std::initializer_list<std::pair<std::initializer_list<size_type>, value_type>> args)
+    /**
+     * Retrieve const reference to the value at the given index.
+     *
+     * @example
+     * mdarray<int> x({1, 1}); // rank 2
+     * x.at({0, 0});           // get the value at (0, 0) in a 2D array.
+     *
+     * @param n Index list
+     * @return reference 
+     */
+    const_reference operator[](const boost::span<const size_type> n) const noexcept
     {
-        for (const auto& arg : args) {
-            this->insert(arg.first, arg.second);
-        }
+        return this->m_data[this->index(n)];
     }
 
     /**
-     * Allocate this mdarray based on the indices given.
+     * Insert a multi-dimensonal value at the given index.
      * 
-     * @example
-     * mdarray<double> x(3);   // rank 3
-     * x.allocate({5, 10, 10}) // Allocates 5 * 10 * 10 elements
-     * @param n 
+     * @param n A multi-dimensional index with the same size as `m_shape`.
+     * @param value The value to set at the given index.
      */
-    void allocate(const boost::span<const size_t> n)
+    void insert(const boost::span<const size_type> n, value_type value)
+    {
+        this->bounds_check(n);
+
+        this->m_data[this->index(n)] = value;
+    }
+
+    void insert(std::initializer_list<std::pair<const boost::span<const size_type>, value_type>> args)
     {
-        this->m_data.resize(this->index(n));
+        for (const auto& arg : args) {
+            this->insert(arg.first, arg.second);
+        }
     }
 
     /**
-     * Get the (total) size of this mdarray
-     * (aka the total number of elements).
-     * 
+     * Get the size of allocated values in this mdarray
+     * (aka the allocated elements in the backing vector)
      * @return size_type
      */
     size_type size() const noexcept
     {
         return this->m_data.size();
     }
 
-    size_type max_size() const noexcept
-    {
-        size_type max = 1;
-        for (const auto& v : this->m_shape) {
-            max *= v;
-        }
-        return max;
-    }
-
     iterator begin() const noexcept {
         return iterator(*this, 0);
     }
@@ -169,13 +177,16 @@ class mdarray
      * @param n 
      * @return size_type 
      */
-    size_type index(const boost::span<const size_t> n) const
+    size_type index(const boost::span<const size_type> n) const
     {
         size_type index = 0, stride = 1;
         for (size_type k = 0; k < this->rank(); k++) {
+            assert(n[k] < m_shape[k]);
             index += n[k] * stride;
             stride *= this->m_shape[k];
         }
+
+        assert(index < this->size());
 
         return index;
     }
@@ -197,6 +208,23 @@ class mdarray
     }
 
   private:
+
+    inline void bounds_check(const boost::span<const size_type> n) const
+    {
+        // This bounds check is required since `index` performs
+        // modulo arithmetic on the dimension indices with carryover,
+        // such that if an mdarray has shape {2, 2}, the index {2, 0}
+        // is equivalent to {0, 1}.
+        for (size_type i = 0; i < this->m_shape.size(); i++) {
+            if (n[i] >= this->m_shape[i])
+                throw std::out_of_range(
+                    "index " + std::to_string(n[i]) +
+                    " must be less than dimension size " +
+                    std::to_string(this->m_shape[i])
+                );
+        }
+    }
+
     std::vector<size_type> m_shape;
     container_type         m_data;
 };

include/utilities/mdframe.hpp

@@ -2,7 +2,5 @@
 #define NGEN_MDFRAME_HPP
 
 #include "mdframe/mdframe.hpp"
-#include "mdframe/handler_csv.hpp"
-#include "mdframe/handler_netcdf.hpp"
 
 #endif // NGEN_MDFRAME_HPP

include/utilities/mdframe/mdframe.hpp

@@ -5,7 +5,7 @@
 #include <unordered_map>
 
 #include "dimension.hpp"
-#include "mdframe/visitors.hpp"
+#include "visitors.hpp"
 #include "variable.hpp"
 
 namespace ngen {

include/utilities/mdframe/variable.hpp

@@ -60,7 +60,7 @@ struct variable {
      */
     variable() noexcept
         : m_name()
-        , m_data(mdarray<int>{0})
+        , m_data()
         , m_dimensions() {};
 
     /**

src/utilities/mdarray/CMakeLists.txt

@@ -0,0 +1,3 @@
+add_library(mdarray INTERFACE)
+add_library(NGen::mdarray ALIAS mdarray)
+target_include_directories(mdarray INTERFACE ${PROJECT_SOURCE_DIR}/include/utilities)

src/utilities/mdframe/CMakeLists.txt

@@ -0,0 +1,4 @@
+add_library(mdframe handler_csv.cpp handler_netcdf.cpp)
+add_library(NGen::mdframe ALIAS mdframe)
+target_include_directories(mdframe PUBLIC ${PROJECT_SOURCE_DIR}/include/utilities)
+target_link_libraries(mdframe PRIVATE NGen::mdarray)

include/utilities/mdframe/handler_csv.hpp → src/utilities/mdframe/handler_csv.cpp

@@ -1,14 +1,44 @@
-#ifndef NGEN_MDFRAME_HANDLERS_CSV_HPP
-#define NGEN_MDFRAME_HANDLERS_CSV_HPP
-
 #include <fstream>
 
-#include "mdframe.hpp"
-#include "cartesian.hpp"
+#include "mdframe/mdframe.hpp"
+
+#include "span.hpp"
 
 namespace ngen {
 
-inline void mdframe::to_csv(const std::string& path, bool header) const
+/**
+ * @brief 
+ * 
+ * @param shape 
+ * @param index 
+ * @param dimension 
+ * @param output 
+ */
+void cartesian_indices_recurse(
+    boost::span<const std::size_t>         shape,
+    std::size_t                            current_dimension_index,
+    std::vector<std::size_t>&              index_buffer,
+    std::vector<std::vector<std::size_t>>& output
+)
+{
+    if (current_dimension_index == shape.size()) {
+        output.push_back(index_buffer);
+        return;
+    }
+
+    for (std::size_t i = 0; i < shape[current_dimension_index]; i++) {
+        index_buffer[current_dimension_index] = i;
+        cartesian_indices_recurse(shape, current_dimension_index + 1, index_buffer, output);
+    }
+}
+
+void cartesian_indices(const boost::span<const std::size_t> shape, std::vector<std::vector<std::size_t>>& output)
+{
+    std::vector<std::size_t> index_buffer(shape.size());
+    cartesian_indices_recurse(shape, 0, index_buffer, output);
+}
+
+void mdframe::to_csv(const std::string& path, bool header) const
 {
     std::ofstream output(path);
     if (!output)
@@ -73,11 +103,10 @@ inline void mdframe::to_csv(const std::string& path, bool header) const
 
     detail::visitors::to_string_visitor visitor{};
     std::string output_line = "";
-    std::vector<size_type> buffer(max_rank, 0);
     std::vector<std::vector<size_type>> indices;
     indices.reserve(rows);
 
-    ngen::cartesian_indices(shape, buffer, 0, indices);
+    ngen::cartesian_indices(shape, indices);
 
     std::vector<size_type> index_buffer(max_rank);
     for (const auto& index : indices) {
@@ -99,4 +128,3 @@ inline void mdframe::to_csv(const std::string& path, bool header) const
 
 } // namespace ngen
 
-#endif // NGEN_MDFRAME_HANDLERS_CSV_HPP