Select product of slices.

include/highfive/bits/H5Slice_traits.hpp

@@ -245,6 +245,71 @@ class HyperSlab {
     std::vector<Select_> selects;
 };
 
+///
+/// \brief Selects the Cartesian product of slices.
+///
+/// Given a one-dimensional dataset one might want to select the union of
+/// multiple, non-overlapping slices. For example,
+///
+///     using Slice = std::array<size_t, 2>;
+///     using Slices = std::vector<Slice>;
+///     auto slices = Slices{{0, 2}, {4, 10}};
+///     dset.select(ProductSet(slices);
+///
+/// to select elements `0`, `1` and `4`, ..., `9` (inclusive).
+///
+/// For a two-dimensional array one which to select the row specified above,
+/// but only columns `2`, `3` and `4`:
+///
+///    dset.select(ProductSet(slices, Slice{2, 5}));
+///    // Analogues with the roles of columns and rows reversed:
+///    dset.select(ProductSet(Slice{2, 5}, slices));
+///
+/// One can generalize once more and allow the unions of slices in both x- and
+/// y-dimension:
+///
+///     auto yslices = Slices{{1, 5}, {7, 8}};
+///     auto xslices = Slices{{0, 3}, {6, 8}};
+///     dset.select(ProductSet(yslices, xslices));
+///
+/// which selects the following from a 11x8 dataset:
+///
+///     . . . . . . . .
+///     x x x . . . x x
+///     x x x . . . x x
+///     x x x . . . x x
+///     x x x . . . x x
+///     . . . . . . . .
+///     . . . . . . . .
+///     x x x . . . x x
+///     . . . . . . . .
+///     . . . . . . . .
+///     . . . . . . . .
+///
+/// Final twist, the selection along one axis may be a point set, from which a
+/// vector of, possibly single-element, slices can be constructed.
+///
+/// The solution works analogous in higher dimensions. A selection `sk` along
+/// axis `k` can be interpreted as a subset `S_k` of the natural numbers. The
+/// index `i` is in `S_k` if it's selected by `sk`.  The `ProductSet` of `s0`,
+/// ..., `sN` selects the Cartesian product `S_0 x ... x S_N`.
+///
+/// Note that the selections along each axis must be sorted and non-overlapping.
+///
+class ProductSet {
+  public:
+    template <class... Slices>
+    explicit ProductSet(const Slices&... slices);
+
+  private:
+    HyperSlab slab;
+    std::vector<size_t> shape;
+
+    template <typename Derivate>
+    friend class SliceTraits;
+};
+
+
 template <typename Derivate>
 class SliceTraits {
   public:
@@ -291,6 +356,8 @@ class SliceTraits {
     ///
     Selection select(const ElementSet& elements) const;
 
+    Selection select(const ProductSet& product_set) const;
+
     template <typename T>
     T read(const DataTransferProps& xfer_props = DataTransferProps()) const;
 

include/highfive/bits/H5Slice_traits_misc.hpp

@@ -63,6 +63,128 @@ inline ElementSet::ElementSet(const std::vector<std::vector<std::size_t>>& eleme
     }
 }
 
+namespace detail {
+class HyperCube {
+  public:
+    HyperCube(size_t rank)
+        : offset(rank)
+        , count(rank) {}
+
+    void cross(const std::array<size_t, 2>& range, size_t axis) {
+        offset[axis] = range[0];
+        count[axis] = range[1] - range[0];
+    }
+
+    RegularHyperSlab asSlab() {
+        return RegularHyperSlab(offset, count);
+    }
+
+  private:
+    std::vector<size_t> offset;
+    std::vector<size_t> count;
+};
+
+inline void build_hyper_slab(HyperSlab& slab, size_t /* axis */, HyperCube& cube) {
+    slab |= cube.asSlab();
+}
+
+template <class... Slices>
+inline void build_hyper_slab(HyperSlab& slab,
+                      size_t axis,
+                      HyperCube& cube,
+                      const std::array<size_t, 2>& slice,
+                      const Slices&... higher_slices) {
+    cube.cross(slice, axis);
+    build_hyper_slab(slab, axis + 1, cube, higher_slices...);
+}
+
+template <class... Slices>
+inline void build_hyper_slab(HyperSlab& slab,
+                      size_t axis,
+                      HyperCube& cube,
+                      const std::vector<std::array<size_t, 2>>& slices,
+                      const Slices&... higher_slices) {
+    for (const auto& slice: slices) {
+        build_hyper_slab(slab, axis, cube, slice, higher_slices...);
+    }
+}
+
+template <class... Slices>
+inline void build_hyper_slab(HyperSlab& slab,
+                      size_t axis,
+                      HyperCube& cube,
+                      const std::vector<size_t>& ids,
+                      const Slices&... higher_slices) {
+    for (const auto& id: ids) {
+        auto slice = std::array<size_t, 2>{id, id + 1};
+        build_hyper_slab(slab, axis, cube, slice, higher_slices...);
+    }
+}
+
+inline void compute_squashed_shape(size_t /* axis */, std::vector<size_t>& /* shape */) {
+    // assert(axis == shape.size());
+}
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                            std::vector<size_t>& shape,
+                            const std::array<size_t, 2>& slice,
+                            const Slices&... higher_slices);
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                            std::vector<size_t>& shape,
+                            const std::vector<size_t>& points,
+                            const Slices&... higher_slices);
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                            std::vector<size_t>& shape,
+                            const std::vector<std::array<size_t, 2>>& slices,
+                            const Slices&... higher_slices);
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                            std::vector<size_t>& shape,
+                            const std::array<size_t, 2>& slice,
+                            const Slices&... higher_slices) {
+    shape[axis] = slice[1] - slice[0];
+    compute_squashed_shape(axis + 1, shape, higher_slices...);
+}
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                            std::vector<size_t>& shape,
+                            const std::vector<size_t>& points,
+                            const Slices&... higher_slices) {
+    shape[axis] = points.size();
+    compute_squashed_shape(axis + 1, shape, higher_slices...);
+}
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                            std::vector<size_t>& shape,
+                            const std::vector<std::array<size_t, 2>>& slices,
+                            const Slices&... higher_slices) {
+    shape[axis] = 0;
+    for (const auto& slice: slices) {
+        shape[axis] += slice[1] - slice[0];
+    }
+    compute_squashed_shape(axis + 1, shape, higher_slices...);
+}
+}  // namespace detail
+
+template <class... Slices>
+inline ProductSet::ProductSet(const Slices&... slices) {
+    auto rank = sizeof...(slices);
+    detail::HyperCube cube(rank);
+    detail::build_hyper_slab(slab, 0, cube, slices...);
+
+    shape = std::vector<size_t>(rank, size_t(0));
+    detail::compute_squashed_shape(0, shape, slices...);
+}
+
+
 template <typename Derivate>
 inline Selection SliceTraits<Derivate>::select(const HyperSlab& hyperslab,
                                                const DataSpace& memspace) const {
@@ -156,6 +278,11 @@ inline Selection SliceTraits<Derivate>::select(const ElementSet& elements) const
     return detail::make_selection(DataSpace(num_elements), space, details::get_dataset(slice));
 }
 
+template <typename Derivate>
+inline Selection SliceTraits<Derivate>::select(const ProductSet& product_set) const {
+    return this->select(product_set.slab, DataSpace(product_set.shape));
+}
+
 
 template <typename Derivate>
 template <typename T>

tests/unit/tests_high_five.hpp

@@ -151,7 +151,7 @@ struct ContentGenerate<std::string> {
         ContentGenerate<char> gen;
         std::string random_string;
         std::mt19937_64 rgen;
-        rgen.seed(88);
+        rgen.seed(42);
         std::uniform_int_distribution<unsigned> int_dist(0, 1000);
         const size_t size_string = int_dist(rgen);
 

tests/unit/tests_high_five_base.cpp

@@ -1617,6 +1617,156 @@ TEMPLATE_LIST_TEST_CASE("irregularHyperSlabSelectionWrite", "[template]", std::t
     irregularHyperSlabSelectionWriteTest<TestType>();
 }
 
+// Ensure that "all" combinations of `ProductSet` are being instantiated.
+class CompileProductSet {
+    using Points = std::vector<size_t>;
+    using Slice = std::array<size_t, 2>;
+    using Slices = std::vector<Slice>;
+
+    void all() {
+        one_arg();
+        two_args();
+        three_args();
+    }
+
+    template <class... Preamble>
+    void zero_args() {
+        ProductSet(Preamble()...);
+    }
+
+    template <class... Preamble>
+    void one_arg() {
+        zero_args<Preamble..., Points>();
+        zero_args<Preamble..., Slice>();
+        zero_args<Preamble..., Slices>();
+    }
+
+    template <class... Preamble>
+    void two_args() {
+        one_arg<Preamble..., Points>();
+        one_arg<Preamble..., Slice>();
+        one_arg<Preamble..., Slices>();
+    }
+
+    template <class... Preamble>
+    void three_args() {
+        two_args<Preamble..., Points>();
+        two_args<Preamble..., Slice>();
+        two_args<Preamble..., Slices>();
+    }
+};
+
+template <class S, class Y, class X>
+void check_product_set_shape(const S& subarray, const Y& yslices, const X& xslices) {
+    std::vector<size_t> subshape{0, 0};
+
+    for (auto yslice: yslices) {
+        subshape[0] += yslice[1] - yslice[0];
+    }
+
+    for (auto xslice: xslices) {
+        subshape[1] += xslice[1] - xslice[0];
+    }
+
+    REQUIRE(subarray.size() == subshape[0]);
+    for (const auto& v: subarray) {
+        REQUIRE(v.size() == subshape[1]);
+    }
+}
+
+template <class S, class Y, class X>
+void check_product_set_values(const S& array,
+                              const S& subarray,
+                              const Y& yslices,
+                              const X& xslices) {
+    size_t il = 0;
+    for (const auto& yslice: yslices) {
+        for (size_t ig = yslice[0]; ig < yslice[1]; ++ig) {
+            size_t jl = 0;
+
+            for (const auto& xslice: xslices) {
+                for (size_t jg = xslice[0]; jg < xslice[1]; ++jg) {
+                    REQUIRE(subarray[il][jl] == array[ig][jg]);
+                    ++jl;
+                }
+            }
+            ++il;
+        }
+    }
+}
+
+template <class S, class Y, class X>
+void check(const S& array, const S& subarray, const Y& yslices, const X& xslices) {
+    check_product_set_shape(subarray, yslices, xslices);
+    check_product_set_values(array, subarray, yslices, xslices);
+};
+
+
+TEST_CASE("productSet") {
+    using Slice = std::array<size_t, 2>;
+    using Slices = std::vector<Slice>;
+    using Points = std::vector<size_t>;
+
+    const std::string file_name("h5_test_product_set.h5");
+
+    auto generate = [](size_t n, size_t m, auto f) {
+        auto x = std::vector<std::vector<double>>(n);
+        for (size_t i = 0; i < n; ++i) {
+            x[i] = std::vector<double>(m);
+        }
+
+        for (size_t i = 0; i < n; ++i) {
+            for (size_t j = 0; j < m; ++j) {
+                x[i][j] = f(i, j);
+            }
+        }
+
+        return x;
+    };
+
+    auto array = generate(6, 12, [](size_t i, size_t j) { return double(i) + double(j) * 0.01; });
+
+    auto file = File(file_name, File::Truncate);
+    auto dset = file.createDataSet("dset", array);
+
+    SECTION("rR") {
+        std::vector<std::vector<double>> subarray;
+
+        auto yslice = Slice{1, 3};
+        auto yslices = Slices{yslice};
+        auto xslices = Slices{{0, 1}, {3, 5}};
+
+        dset.select(ProductSet(yslice, xslices)).read(subarray);
+
+        check(array, subarray, yslices, xslices);
+    }
+
+    SECTION("Rr") {
+        std::vector<std::vector<double>> subarray;
+
+        auto yslices = Slices{{0, 1}, {3, 5}};
+        auto xslice = Slice{1, 3};
+        auto xslices = Slices{xslice};
+
+        dset.select(ProductSet(yslices, xslice)).read(subarray);
+
+        check(array, subarray, yslices, xslices);
+    }
+
+    SECTION("Rp") {
+        std::vector<std::vector<double>> subarray;
+
+        auto yslices = Slices{{0, 1}, {3, 5}};
+        auto xpoints = Points{2, 4, 5};
+        auto xslices = Slices{{2, 3}, {4, 6}};
+
+        dset.select(ProductSet(yslices, xpoints)).read(subarray);
+
+        check(array, subarray, yslices, xslices);
+    }
+}
+
+
 template <typename T>
 void attribute_scalar_rw() {
     std::ostringstream filename;