highfive/bits/H5Slice_traits_misc.hpp

0001 /*
0002  *  Copyright (c), 2017, Adrien Devresse <adrien.devresse@epfl.ch>
0003  *
0004  *  Distributed under the Boost Software License, Version 1.0.
0005  *    (See accompanying file LICENSE_1_0.txt or copy at
0006  *          http://www.boost.org/LICENSE_1_0.txt)
0007  *
0008  */
0009 #pragma once
0010
0011 #include <algorithm>
0012 #include <cassert>
0013 #include <functional>
0014 #include <numeric>
0015 #include <sstream>
0016 #include <string>
0017
0018 #include "h5d_wrapper.hpp"
0019 #include "h5s_wrapper.hpp"
0020
0021 #include "H5ReadWrite_misc.hpp"
0022 #include "H5Converter_misc.hpp"
0023 #include "squeeze.hpp"
0024 #include "compute_total_size.hpp"
0025 #include "assert_compatible_spaces.hpp"
0026
0027 namespace HighFive {
0028
0029 namespace details {
0030
0031 // map the correct reference to the dataset depending of the layout
0032 // dataset -> itself
0033 // subselection -> parent dataset
0034 inline const DataSet& get_dataset(const Selection& sel) {
0035     return sel.getDataset();
0036 }
0037
0038 inline const DataSet& get_dataset(const DataSet& ds) {
0039     return ds;
0040 }
0041
0042 // map the correct memspace identifier depending of the layout
0043 // dataset -> entire memspace
0044 // selection -> resolve space id
0045 inline hid_t get_memspace_id(const Selection& ptr) {
0046     return ptr.getMemSpace().getId();
0047 }
0048
0049 inline hid_t get_memspace_id(const DataSet&) {
0050     return H5S_ALL;
0051 }
0052 }  // namespace details
0053
0054 inline ElementSet::ElementSet(std::initializer_list<std::size_t> list)
0055     : _ids(list) {}
0056
0057 inline ElementSet::ElementSet(std::initializer_list<std::vector<std::size_t>> list)
0058     : ElementSet(std::vector<std::vector<std::size_t>>(list)) {}
0059
0060 inline ElementSet::ElementSet(const std::vector<std::size_t>& element_ids)
0061     : _ids(element_ids) {}
0062
0063 inline ElementSet::ElementSet(const std::vector<std::vector<std::size_t>>& element_ids) {
0064     for (const auto& vec: element_ids) {
0065         std::copy(vec.begin(), vec.end(), std::back_inserter(_ids));
0066     }
0067 }
0068
0069 namespace detail {
0070 class HyperCube {
0071   public:
0072     explicit HyperCube(size_t rank)
0073         : offset(rank)
0074         , count(rank) {}
0075
0076     void cross(const std::array<size_t, 2>& range, size_t axis) {
0077         offset[axis] = range[0];
0078         count[axis] = range[1] - range[0];
0079     }
0080
0081     RegularHyperSlab asSlab() {
0082         return RegularHyperSlab(offset, count);
0083     }
0084
0085   private:
0086     std::vector<size_t> offset;
0087     std::vector<size_t> count;
0088 };
0089
0090 inline void build_hyper_slab(HyperSlab& slab, size_t /* axis */, HyperCube& cube) {
0091     slab |= cube.asSlab();
0092 }
0093
0094 template <class... Slices>
0095 inline void build_hyper_slab(HyperSlab& slab,
0096                              size_t axis,
0097                              HyperCube& cube,
0098                              const std::array<size_t, 2>& slice,
0099                              const Slices&... higher_slices) {
0100     cube.cross(slice, axis);
0101     build_hyper_slab(slab, axis + 1, cube, higher_slices...);
0102 }
0103
0104 template <class... Slices>
0105 inline void build_hyper_slab(HyperSlab& slab,
0106                              size_t axis,
0107                              HyperCube& cube,
0108                              const std::vector<std::array<size_t, 2>>& slices,
0109                              const Slices&... higher_slices) {
0110     for (const auto& slice: slices) {
0111         build_hyper_slab(slab, axis, cube, slice, higher_slices...);
0112     }
0113 }
0114
0115 template <class... Slices>
0116 inline void build_hyper_slab(HyperSlab& slab,
0117                              size_t axis,
0118                              HyperCube& cube,
0119                              const std::vector<size_t>& ids,
0120                              const Slices&... higher_slices) {
0121     for (const auto& id: ids) {
0122         auto slice = std::array<size_t, 2>{id, id + 1};
0123         build_hyper_slab(slab, axis, cube, slice, higher_slices...);
0124     }
0125 }
0126
0127 template <class... Slices>
0128 inline void build_hyper_slab(HyperSlab& slab,
0129                              size_t axis,
0130                              HyperCube& cube,
0131                              size_t id,
0132                              const Slices&... higher_slices) {
0133     auto slice = std::array<size_t, 2>{id, id + 1};
0134     build_hyper_slab(slab, axis, cube, slice, higher_slices...);
0135 }
0136
0137 inline void compute_squashed_shape(size_t /* axis */, std::vector<size_t>& /* shape */) {
0138     // assert(axis == shape.size());
0139 }
0140
0141 template <class... Slices>
0142 inline void compute_squashed_shape(size_t axis,
0143                                    std::vector<size_t>& shape,
0144                                    const std::array<size_t, 2>& slice,
0145                                    const Slices&... higher_slices);
0146
0147 template <class... Slices>
0148 inline void compute_squashed_shape(size_t axis,
0149                                    std::vector<size_t>& shape,
0150                                    const std::vector<size_t>& points,
0151                                    const Slices&... higher_slices);
0152
0153 template <class... Slices>
0154 inline void compute_squashed_shape(size_t axis,
0155                                    std::vector<size_t>& shape,
0156                                    size_t point,
0157                                    const Slices&... higher_slices);
0158
0159 template <class... Slices>
0160 inline void compute_squashed_shape(size_t axis,
0161                                    std::vector<size_t>& shape,
0162                                    const std::vector<std::array<size_t, 2>>& slices,
0163                                    const Slices&... higher_slices);
0164
0165 template <class... Slices>
0166 inline void compute_squashed_shape(size_t axis,
0167                                    std::vector<size_t>& shape,
0168                                    const std::array<size_t, 2>& slice,
0169                                    const Slices&... higher_slices) {
0170     shape[axis] = slice[1] - slice[0];
0171     compute_squashed_shape(axis + 1, shape, higher_slices...);
0172 }
0173
0174 template <class... Slices>
0175 inline void compute_squashed_shape(size_t axis,
0176                                    std::vector<size_t>& shape,
0177                                    const std::vector<size_t>& points,
0178                                    const Slices&... higher_slices) {
0179     shape[axis] = points.size();
0180     compute_squashed_shape(axis + 1, shape, higher_slices...);
0181 }
0182
0183 template <class... Slices>
0184 inline void compute_squashed_shape(size_t axis,
0185                                    std::vector<size_t>& shape,
0186                                    const std::vector<std::array<size_t, 2>>& slices,
0187                                    const Slices&... higher_slices) {
0188     shape[axis] = 0;
0189     for (const auto& slice: slices) {
0190         shape[axis] += slice[1] - slice[0];
0191     }
0192     compute_squashed_shape(axis + 1, shape, higher_slices...);
0193 }
0194
0195 template <class... Slices>
0196 inline void compute_squashed_shape(size_t axis,
0197                                    std::vector<size_t>& shape,
0198                                    size_t /* point */,
0199                                    const Slices&... higher_slices) {
0200     shape[axis] = 1;
0201     compute_squashed_shape(axis + 1, shape, higher_slices...);
0202 }
0203 }  // namespace detail
0204
0205 template <class... Slices>
0206 inline ProductSet::ProductSet(const Slices&... slices) {
0207     auto rank = sizeof...(slices);
0208     detail::HyperCube cube(rank);
0209     detail::build_hyper_slab(slab, 0, cube, slices...);
0210
0211     shape = std::vector<size_t>(rank, size_t(0));
0212     detail::compute_squashed_shape(0, shape, slices...);
0213 }
0214
0215
0216 template <typename Derivate>
0217 inline Selection SliceTraits<Derivate>::select(const HyperSlab& hyper_slab,
0218                                                const DataSpace& memspace) const {
0219     // Note: The current limitation are that memspace must describe a
0220     //       packed memspace.
0221     //
0222     //       The reason for this is that we're unable to unpack general
0223     //       hyperslabs when the memory is not contiguous, e.g.
0224     //       `std::vector<std::vector<double>>`.
0225     const auto& slice = static_cast<const Derivate&>(*this);
0226     auto filespace = hyper_slab.apply(slice.getSpace());
0227
0228     return detail::make_selection(memspace, filespace, details::get_dataset(slice));
0229 }
0230
0231 template <typename Derivate>
0232 template <typename Impl>
0233 inline Selection SliceTraits<Derivate>::select(const HyperSlabInterface<Impl>& hyper_slab) const {
0234     const auto& slice = static_cast<const Derivate&>(*this);
0235     auto filespace = slice.getSpace();
0236     filespace = hyper_slab.apply(filespace);
0237
0238     const auto n_elements = detail::h5s_get_select_npoints(filespace.getId());
0239     auto memspace = DataSpace{static_cast<size_t>(n_elements)};
0240
0241     return detail::make_selection(memspace, filespace, details::get_dataset(slice));
0242 }
0243
0244 template <typename Derivate>
0245 inline Selection SliceTraits<Derivate>::select(const std::vector<size_t>& offset,
0246                                                const std::vector<size_t>& count,
0247                                                const std::vector<size_t>& stride,
0248                                                const std::vector<size_t>& block) const {
0249     auto slab = HyperSlab(RegularHyperSlab(offset, count, stride, block));
0250     auto memspace = DataSpace(count);
0251     return select(slab, memspace);
0252 }
0253
0254 template <typename Derivate>
0255 inline Selection SliceTraits<Derivate>::select(const std::vector<size_t>& columns) const {
0256     const auto& slice = static_cast<const Derivate&>(*this);
0257     const DataSpace& space = slice.getSpace();
0258     std::vector<size_t> dims = space.getDimensions();
0259
0260     if (dims.empty()) {
0261         throw DataSpaceException(
0262             "Invalid, zero-dimensional (scalar) dataspace encountered when "
0263             "selecting columns; must be atleast 1-dimensional.");
0264     }
0265
0266     std::vector<size_t> counts = dims;
0267     counts.back() = 1;
0268
0269     std::vector<size_t> offsets(dims.size(), 0);
0270
0271     HyperSlab slab;
0272     for (const auto& column: columns) {
0273         offsets.back() = column;
0274         slab |= RegularHyperSlab(offsets, counts);
0275     }
0276
0277     std::vector<size_t> memdims = dims;
0278     memdims.back() = columns.size();
0279
0280     return select(slab, DataSpace(memdims));
0281 }
0282
0283 template <typename Derivate>
0284 inline Selection SliceTraits<Derivate>::select(const ElementSet& elements) const {
0285     const auto& slice = static_cast<const Derivate&>(*this);
0286     const hsize_t* data = nullptr;
0287     const DataSpace space = slice.getSpace().clone();
0288     const std::size_t length = elements._ids.size();
0289     if (length % space.getNumberDimensions() != 0) {
0290         throw DataSpaceException(
0291             "Number of coordinates in elements picking "
0292             "should be a multiple of the dimensions.");
0293     }
0294     const std::size_t num_elements = length / space.getNumberDimensions();
0295     std::vector<hsize_t> raw_elements;
0296
0297     // optimised at compile time
0298     // switch for data conversion on 32bits platforms
0299     if (std::is_same<std::size_t, hsize_t>::value) {
0300         // `if constexpr` can't be used, thus a reinterpret_cast is needed.
0301         data = reinterpret_cast<const hsize_t*>(elements._ids.data());
0302     } else {
0303         raw_elements.resize(length);
0304         std::copy(elements._ids.begin(), elements._ids.end(), raw_elements.begin());
0305         data = raw_elements.data();
0306     }
0307
0308     detail::h5s_select_elements(space.getId(), H5S_SELECT_SET, num_elements, data);
0309
0310     return detail::make_selection(DataSpace(num_elements), space, details::get_dataset(slice));
0311 }
0312
0313 template <typename Derivate>
0314 inline Selection SliceTraits<Derivate>::select(const ProductSet& product_set) const {
0315     return this->select(product_set.slab, DataSpace(product_set.shape));
0316 }
0317
0318
0319 template <typename Derivate>
0320 template <typename T>
0321 inline T SliceTraits<Derivate>::read(const DataTransferProps& xfer_props) const {
0322     T array;
0323     read(array, xfer_props);
0324     return array;
0325 }
0326
0327
0328 template <typename Derivate>
0329 template <typename T>
0330 inline void SliceTraits<Derivate>::read(T& array, const DataTransferProps& xfer_props) const {
0331     const auto& slice = static_cast<const Derivate&>(*this);
0332     const DataSpace& mem_space = slice.getMemSpace();
0333
0334     auto file_datatype = slice.getDataType();
0335
0336     const details::BufferInfo<T> buffer_info(
0337         file_datatype,
0338         [&slice]() -> std::string { return details::get_dataset(slice).getPath(); },
0339         details::BufferInfo<T>::Operation::read);
0340
0341     if (!details::checkDimensions(mem_space, buffer_info.getMinRank(), buffer_info.getMaxRank())) {
0342         std::ostringstream ss;
0343         ss << "Impossible to read DataSet of dimensions " << mem_space.getNumberDimensions()
0344            << " into arrays of dimensions: " << buffer_info.getMinRank() << "(min) to "
0345            << buffer_info.getMaxRank() << "(max)";
0346         throw DataSpaceException(ss.str());
0347     }
0348     auto dims = mem_space.getDimensions();
0349
0350     auto r = details::data_converter::get_reader<T>(dims, array, file_datatype);
0351     read_raw(r.getPointer(), buffer_info.data_type, xfer_props);
0352     // re-arrange results
0353     r.unserialize(array);
0354
0355     auto t = buffer_info.data_type;
0356     auto c = t.getClass();
0357     if (c == DataTypeClass::VarLen || t.isVariableStr()) {
0358 #if H5_VERSION_GE(1, 12, 0)
0359         // This one have been created in 1.12.0
0360         (void)
0361             detail::h5t_reclaim(t.getId(), mem_space.getId(), xfer_props.getId(), r.getPointer());
0362 #else
0363         // This one is deprecated since 1.12.0
0364         (void) detail::h5d_vlen_reclaim(t.getId(),
0365                                         mem_space.getId(),
0366                                         xfer_props.getId(),
0367                                         r.getPointer());
0368 #endif
0369     }
0370 }
0371
0372
0373 template <typename Derivate>
0374 template <typename T>
0375 inline void SliceTraits<Derivate>::read_raw(T* array,
0376                                             const DataType& mem_datatype,
0377                                             const DataTransferProps& xfer_props) const {
0378     static_assert(!std::is_const<T>::value,
0379                   "read() requires a non-const structure to read data into");
0380
0381     const auto& slice = static_cast<const Derivate&>(*this);
0382
0383     detail::h5d_read(details::get_dataset(slice).getId(),
0384                      mem_datatype.getId(),
0385                      details::get_memspace_id(slice),
0386                      slice.getSpace().getId(),
0387                      xfer_props.getId(),
0388                      static_cast<void*>(array));
0389 }
0390
0391
0392 template <typename Derivate>
0393 template <typename T>
0394 inline void SliceTraits<Derivate>::read_raw(T* array, const DataTransferProps& xfer_props) const {
0395     using element_type = typename details::inspector<T>::base_type;
0396     const DataType& mem_datatype = create_and_check_datatype<element_type>();
0397
0398     read_raw(array, mem_datatype, xfer_props);
0399 }
0400
0401
0402 template <typename Derivate>
0403 template <typename T>
0404 inline void SliceTraits<Derivate>::write(const T& buffer, const DataTransferProps& xfer_props) {
0405     const auto& slice = static_cast<const Derivate&>(*this);
0406     const DataSpace& mem_space = slice.getMemSpace();
0407     auto dims = mem_space.getDimensions();
0408
0409     auto file_datatype = slice.getDataType();
0410
0411     const details::BufferInfo<T> buffer_info(
0412         file_datatype,
0413         [&slice]() -> std::string { return details::get_dataset(slice).getPath(); },
0414         details::BufferInfo<T>::Operation::write);
0415
0416     if (!details::checkDimensions(mem_space, buffer_info.getMinRank(), buffer_info.getMaxRank())) {
0417         std::ostringstream ss;
0418         ss << "Impossible to write buffer with dimensions n = " << buffer_info.getRank(buffer)
0419            << "into dataset with dimensions " << details::format_vector(mem_space.getDimensions())
0420            << ".";
0421         throw DataSpaceException(ss.str());
0422     }
0423     auto w = details::data_converter::serialize<T>(buffer, dims, file_datatype);
0424     write_raw(w.getPointer(), buffer_info.data_type, xfer_props);
0425 }
0426
0427
0428 template <typename Derivate>
0429 template <typename T>
0430 inline void SliceTraits<Derivate>::write_raw(const T* buffer,
0431                                              const DataType& mem_datatype,
0432                                              const DataTransferProps& xfer_props) {
0433     const auto& slice = static_cast<const Derivate&>(*this);
0434
0435     detail::h5d_write(details::get_dataset(slice).getId(),
0436                       mem_datatype.getId(),
0437                       details::get_memspace_id(slice),
0438                       slice.getSpace().getId(),
0439                       xfer_props.getId(),
0440                       static_cast<const void*>(buffer));
0441 }
0442
0443
0444 template <typename Derivate>
0445 template <typename T>
0446 inline void SliceTraits<Derivate>::write_raw(const T* buffer, const DataTransferProps& xfer_props) {
0447     using element_type = typename details::inspector<T>::base_type;
0448     const auto& mem_datatype = create_and_check_datatype<element_type>();
0449
0450     write_raw(buffer, mem_datatype, xfer_props);
0451 }
0452
0453 namespace detail {
0454 inline const DataSet& getDataSet(const Selection& selection) {
0455     return selection.getDataset();
0456 }
0457
0458 inline const DataSet& getDataSet(const DataSet& dataset) {
0459     return dataset;
0460 }
0461
0462 }  // namespace detail
0463
0464 template <typename Derivate>
0465 inline Selection SliceTraits<Derivate>::squeezeMemSpace(const std::vector<size_t>& axes) const {
0466     auto slice = static_cast<const Derivate&>(*this);
0467     auto mem_dims = slice.getMemSpace().getDimensions();
0468     auto squeezed_dims = detail::squeeze(mem_dims, axes);
0469
0470     return detail::make_selection(DataSpace(squeezed_dims),
0471                                   slice.getSpace(),
0472                                   detail::getDataSet(slice));
0473 }
0474
0475 template <typename Derivate>
0476 inline Selection SliceTraits<Derivate>::reshapeMemSpace(const std::vector<size_t>& new_dims) const {
0477     auto slice = static_cast<const Derivate&>(*this);
0478
0479     detail::assert_compatible_spaces(slice.getMemSpace(), new_dims);
0480     return detail::make_selection(DataSpace(new_dims), slice.getSpace(), detail::getDataSet(slice));
0481 }
0482
0483 }  // namespace HighFive