33#error "HDF5 support is not enabled in this build of nda. Please configure and install nda with -DHDF5Support=ON"
46 template <MemoryArray A>
47 void resize_or_check(A &a, std::array<long, A::rank>
const &shape) {
48 if constexpr (is_regular_v<A>) {
51 if (a.shape() != shape) NDA_RUNTIME_ERROR <<
"Error in nda::detail::resize_or_check: Dimension mismatch: " << shape <<
" != " << a.shape();
56 template <MemoryArray A>
57 auto prepare_h5_array_view(
const A &a) {
58 auto [parent_shape, h5_strides] =
59 h5::array_interface::get_parent_shape_and_h5_strides(a.indexmap().strides().data(), A::rank, a.shape().data());
61 for (
int u = 0; u < A::rank; ++u) {
62 v.slab.count[u] = a.shape()[u];
63 v.slab.stride[u] = h5_strides[u];
64 v.parent_shape[u] = parent_shape[u];
70 template <MemoryArrayOfRank<1> A>
71 h5::char_buf to_char_buf(A
const &a) {
74 for (
auto &x : a) s = std::max(s, x.size() + 1);
77 std::vector<char> buf;
78 buf.resize(a.size() * s, 0x00);
81 strcpy(&buf[i * s], x.c_str());
86 auto len = h5::v_t{size_t(a.size()), s};
91 template <MemoryArrayOfRank<1> A>
92 void from_char_buf(h5::char_buf
const &cb, A &a) {
94 resize_or_check(a, std::array<long, 1>{
static_cast<long>(cb.lengths[0])});
97 auto len_string = cb.lengths[1];
101 x.append(&cb.buffer[i * len_string]);
123 template <MemoryArray A>
124 void h5_write(h5::group g, std::string
const &name, A
const &a,
bool compress =
true) {
125 if constexpr (std::is_same_v<nda::get_value_t<A>, std::string>) {
127 h5_write(g, name, detail::to_char_buf(a));
131 if (not a.indexmap().is_stride_order_C()) {
133 auto a_c_layout = h5_arr_t{a.shape()};
135 h5_write(g, name, a_c_layout, compress);
140 auto v = detail::prepare_h5_array_view(a);
141 h5::array_interface::write(g, name, v, compress);
144 auto g2 = g.create_group(name);
147 nda::for_each(a.shape(), [&](
auto... is) { h5_write(g2, make_name(is...), a(is...)); });
166 template <
size_t NDim,
typename... IRs>
169 static constexpr auto size_of_slice =
sizeof...(IRs);
172 static constexpr auto ellipsis_count = (std::is_same_v<IRs, ellipsis> + ... + 0);
173 static_assert(ellipsis_count < 2,
"Error in nda::hyperslab_and_shape_from_slice: Only a single ellipsis is allowed in slicing");
174 static constexpr auto has_ellipsis = (ellipsis_count == 1);
177 static constexpr auto ellipsis_position = [&]<
size_t... Is>(std::index_sequence<Is...>) {
178 if constexpr (has_ellipsis)
return ((std::is_same_v<IRs, ellipsis> * Is) + ... + 0);
179 return size_of_slice;
180 }(std::index_sequence_for<IRs...>{});
183 static constexpr auto integer_count = (std::integral<IRs> + ... + 0);
186 static constexpr auto range_count = size_of_slice - integer_count - ellipsis_count;
187 static_assert((has_ellipsis && range_count <= NDim) || range_count == NDim,
188 "Error in nda::hyperslab_and_shape_from_slice: Rank does not match the number of non-trivial slice dimensions");
191 static constexpr auto ellipsis_width = NDim - range_count;
194 static constexpr auto slab_rank = NDim + integer_count;
197 auto ds_rank = ds_shape.size() - is_complex;
198 if (slab_rank != ds_rank)
199 NDA_RUNTIME_ERROR <<
"Error in nda::hyperslab_and_shape_from_slice: Incompatible dataset and slice ranks: " << ds_rank
200 <<
" != " << size_of_slice;
203 auto slab = h5::array_interface::hyperslab(slab_rank, is_complex);
204 auto shape = std::array<long, NDim>{};
205 [&, m = 0]<
size_t... Is>(std::index_sequence<Is...>)
mutable {
207 [&]<
typename IR>(
size_t n, IR
const &ir)
mutable {
208 if (n > ellipsis_position) n += (ellipsis_width - 1);
209 if constexpr (std::integral<IR>) {
212 }
else if constexpr (std::is_same_v<IR, nda::ellipsis>) {
213 for (
auto k : range(n, n + ellipsis_width)) {
214 slab.count[k] = ds_shape[k];
215 shape[m++] = ds_shape[k];
217 }
else if constexpr (std::is_same_v<IR, nda::range>) {
218 slab.offset[n] = ir.first();
219 slab.stride[n] = ir.step();
220 slab.count[n] = ir.size();
221 shape[m++] = ir.size();
223 static_assert(std::is_same_v<IR, nda::range::all_t>);
224 slab.count[n] = ds_shape[n];
225 shape[m++] = ds_shape[n];
227 }(Is, std::get<Is>(slice)),
229 }(std::make_index_sequence<size_of_slice>{});
230 return std::make_pair(slab, shape);
249 template <MemoryArray A,
typename... IRs>
250 void h5_write(h5::group g, std::string
const &name, A
const &a, std::tuple<IRs...>
const &slice) {
252 constexpr int size_of_slice =
sizeof...(IRs);
253 static_assert(size_of_slice > 0,
"Error in nda::h5_write: Invalid empty slice");
258 if (not a.indexmap().is_stride_order_C()) {
260 auto a_c_layout = h5_arr_t{a.shape()};
262 h5_write(g, name, a_c_layout, slice);
267 auto ds_info = h5::array_interface::get_dataset_info(g, name);
268 if (is_complex != ds_info.has_complex_attribute)
269 NDA_RUNTIME_ERROR <<
"Error in nda::h5_write: Dataset and array/view must both be complex or non-complex";
273 if (sh != a.shape()) NDA_RUNTIME_ERROR <<
"Error in nda::h5_write: Incompatible slice and array shape: " << sh <<
" != " << a.shape();
276 auto v = detail::prepare_h5_array_view(a);
277 h5::array_interface::write_slice(g, name, v, sl);
301 template <MemoryArray A,
typename... IRs>
302 void h5_read(h5::group g, std::string
const &name, A &a, std::tuple<IRs...>
const &slice = {}) {
303 constexpr bool slicing = (
sizeof...(IRs) > 0);
305 if constexpr (std::is_same_v<typename A::value_type, std::string>) {
307 static_assert(!slicing,
"Error in nda::h5_read: Slicing not supported for arrays/views of strings");
310 detail::from_char_buf(cb, a);
314 if (not a.indexmap().is_stride_order_C()) {
316 auto a_c_layout = h5_arr_t{};
317 h5_read(g, name, a_c_layout, slice);
318 detail::resize_or_check(a, a_c_layout.shape());
324 auto ds_info = h5::array_interface::get_dataset_info(g, name);
328 if constexpr (is_complex) {
329 if (!ds_info.has_complex_attribute) {
338 std::array<long, A::rank> shape{};
339 auto slab = h5::array_interface::hyperslab{};
340 if constexpr (slicing) {
345 for (
int u = 0; u < A::rank; ++u) shape[u] = ds_info.lengths[u];
349 detail::resize_or_check(a, shape);
352 auto ds_rank = ds_info.rank() - is_complex;
353 if (!slicing && ds_rank != A::rank)
354 NDA_RUNTIME_ERROR <<
"Error in nda::h5_read: Incompatible dataset and array ranks: " << ds_rank <<
" != " << A::rank;
357 auto v = detail::prepare_h5_array_view(a);
358 h5::array_interface::read(g, name, v, slab);
361 static_assert(!slicing,
"Error in nda::h5_read: Slicing not supported for arrays/views of generic types");
362 auto g2 = g.open_group(name);
365 std::array<long, A::rank> h5_shape;
366 h5_read(g2,
"shape", h5_shape);
367 detail::resize_or_check(a, h5_shape);
371 nda::for_each(a.shape(), [&](
auto... is) { h5_read(g2, make_name(is...), a(is...)); });
Provides concepts for the nda library.
Provides various convenient aliases and helper functions for nda::basic_array and nda::basic_array_vi...
Provides a custom runtime error class and macros to assert conditions and throw exceptions.
Provides for_each functions for multi-dimensional arrays/views.
void h5_read(h5::group g, std::string const &name, A &a, std::tuple< IRs... > const &slice={})
Read into an nda::MemoryArray from an HDF5 file.
void h5_write(h5::group g, std::string const &name, A const &a, bool compress=true)
Write an nda::MemoryArray to a new dataset/subgroup into an HDF5 file.
auto hyperslab_and_shape_from_slice(std::tuple< IRs... > const &slice, std::vector< h5::hsize_t > const &ds_shape, bool is_complex)
Construct an h5::array_interface::hyperslab and the corresponding shape from a given slice,...
basic_array< ValueType, Rank, Layout, 'A', ContainerPolicy > array
Alias template of an nda::basic_array with an 'A' algebra.
__inline__ void for_each(std::array< Int, R > const &shape, F &&f)
Loop over all possible index values of a given shape and apply a function to them.
std::string to_string(std::array< T, R > const &a)
Get a string representation of a std::array.
constexpr bool is_complex_v
Constexpr variable that is true if type T is a std::complex type.
constexpr bool is_scalar_v
Constexpr variable that is true if type S is a scalar type, i.e. arithmetic or complex.
Includes the itertools header and provides some additional utilities.
Provides type traits for the nda library.