Select product of slices.

include/highfive/bits/H5Slice_traits.hpp

@@ -231,6 +231,85 @@ 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 and axis may be discrete indices, from
+/// which a vector of, possibly single-element, slices can be constructed. The
+/// corresponding types are `std::vector<size_t>` and `size_t` for multiple or
+/// just a single values. Note, looping over rows or columns one-by-one can be
+/// a very serious performance problem. In particular,
+///
+///     // Avoid:
+///     for(auto i : indices) {
+///         dset.select(i).read<double>();
+///     }
+///
+///     // Use:
+///     std::vector<size_t> tmp(indices.begin(), indices.end());
+///     dset.select(tmp).read<std::vector<double>>();
+///
+/// obvious omit the copy if `indices` already has the correct type.
+///
+/// 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:
@@ -277,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

@@ -66,6 +66,153 @@ 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...);
+    }
+}
+
+template <class... Slices>
+inline void build_hyper_slab(HyperSlab& slab,
+                             size_t axis,
+                             HyperCube& cube,
+                             size_t id,
+                             const Slices&... higher_slices) {
+    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,
+                                   size_t point,
+                                   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...);
+}
+
+template <class... Slices>
+inline void compute_squashed_shape(size_t axis,
+                                   std::vector<size_t>& shape,
+                                   size_t /* point */,
+                                   const Slices&... higher_slices) {
+    shape[axis] = 1;
+    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 {
@@ -157,6 +304,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>

src/examples/select_slices.cpp

@@ -0,0 +1,131 @@
+#include <functional>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include <highfive/highfive.hpp>
+
+void print_mask(const std::vector<std::vector<double>>& values,
+                const std::vector<std::vector<double>>& selected);
+
+void print_values(const std::vector<std::vector<double>>& values);
+
+void pretty_print(const std::vector<std::vector<double>>& values,
+                  const std::vector<std::vector<double>>& selected);
+
+int main(void) {
+    using namespace HighFive;
+    using container_type = std::vector<std::vector<double>>;
+
+    const std::string file_name("select_slices.h5");
+    const std::string dataset_name("dset");
+
+    // Create a new file using the default property lists.
+    File file(file_name, File::Truncate);
+
+    // we have some example values in a 2D vector 2x5
+    // clang-format off
+    container_type values = {
+      {1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7},
+      {2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7},
+      {3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7},
+      {4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7},
+      {5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7},
+      {6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7},
+      {7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7},
+      {8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7},
+      {9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7}
+    };
+    // clang-format on
+
+    auto dset = file.createDataSet(dataset_name, values);
+
+    // Let's start with the selection `values[1:4, 2:4]`.
+    {
+        auto xslice = std::array<size_t, 2>{2, 4};
+        auto yslice = std::array<size_t, 2>{1, 4};
+
+        auto selected = dset.select(ProductSet(yslice, xslice)).read<container_type>();
+        std::cout << " -- values[1:4, 2:4] ------------ \n";
+        pretty_print(values, selected);
+    }
+
+    // Now we'd like the selection `values[[1,2,8], 2:4]`.
+    {
+        auto xslice = std::array<size_t, 2>{2, 4};
+        auto yslice = std::vector<size_t>{1, 2, 8};
+
+        auto selected = dset.select(ProductSet(yslice, xslice)).read<container_type>();
+        std::cout << "\n -- values[[1,2,8], 2:4] -------- \n";
+        pretty_print(values, selected);
+    }
+
+    // ... or the union of multiple slices.
+    {
+        auto xslice = std::array<size_t, 2>{2, 4};
+        auto yslice = std::vector<std::array<size_t, 2>>{{0, 2}, {5, 9}};
+
+        auto selected = dset.select(ProductSet(yslice, xslice)).read<container_type>();
+        std::cout << "\n -- values[[0:2, 5:10], 2:4] -------- \n";
+        pretty_print(values, selected);
+    }
+
+    // Similar for the union of multiple slices in both x- and y-direction.
+    {
+        auto xslice = std::vector<std::array<size_t, 2>>{{0, 1}, {2, 4}, {6, 7}};
+        auto yslice = std::vector<std::array<size_t, 2>>{{0, 2}, {5, 9}};
+
+        auto selected = dset.select(ProductSet(yslice, xslice)).read<container_type>();
+        std::cout << "\n -- values[[0:2, 5:10], [0:1, 2:4, 6:7]] -------- \n";
+        pretty_print(values, selected);
+    }
+
+    // If only a single row/column is need use the following. However,
+    // selecting elements one-by-one in a loop can be a serious performance
+    // issue.
+    {
+        auto xslice = std::vector<std::array<size_t, 2>>{{0, 1}, {2, 4}, {6, 7}};
+        auto row_id = 3;
+
+        auto selected = dset.select(ProductSet(row_id, xslice)).read<container_type>();
+        std::cout << "\n -- values[3, [0:1, 2:4, 6:7]] -------- \n";
+        pretty_print(values, selected);
+    }
+
+    return 0;
+}
+
+void print_mask(const std::vector<std::vector<double>>& values,
+                const std::vector<std::vector<double>>& selected) {
+    std::vector<double> flat_selection;
+    for (const auto& row: selected) {
+        for (const auto& x: row) {
+            flat_selection.push_back(x);
+        }
+    }
+
+    for (const auto& row: values) {
+        for (const auto& x: row) {
+            auto found = std::find(flat_selection.begin(), flat_selection.end(), x) !=
+                         flat_selection.end();
+            std::cout << (found ? "x" : ".") << "  ";
+        }
+        std::cout << "\n";
+    }
+}
+
+void print_values(const std::vector<std::vector<double>>& values) {
+    for (const auto& row: values) {
+        for (const auto& x: row) {
+            std::cout << x << "  ";
+        }
+        std::cout << "\n";
+    }
+}
+
+void pretty_print(const std::vector<std::vector<double>>& values,
+                  const std::vector<std::vector<double>>& selected) {
+    print_mask(values, selected);
+    std::cout << "\n";
+    print_values(selected);
+}