TRIQS/nda 1.3.0
Multi-dimensional array library for C++
Loading...
Searching...
No Matches
h5.hpp
Go to the documentation of this file.
1// Copyright (c) 2019--present, The Simons Foundation
2// This file is part of TRIQS/nda and is licensed under the Apache License, Version 2.0.
3// SPDX-License-Identifier: Apache-2.0
4// See LICENSE in the root of this distribution for details.
5
10
11#pragma once
12
13#include "./concepts.hpp"
14#include "./declarations.hpp"
15#include "./exceptions.hpp"
16#include "./layout/for_each.hpp"
17#include "./layout/range.hpp"
18#include "./traits.hpp"
19
20#include <h5/h5.hpp>
21
22#include <algorithm>
23#include <array>
24#include <concepts>
25#include <cstddef>
26#include <string>
27#include <tuple>
28#include <type_traits>
29#include <utility>
30#include <vector>
31
32#ifndef NDA_HAVE_H5
33#error "HDF5 support is not enabled in this build of nda. Please configure and install nda with -DHDF5Support=ON"
34#endif
35
36namespace nda {
37
42
43 namespace detail {
44
45 // Resize a given array to fit a given shape or check if a given view fits the shape.
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>) {
49 a.resize(shape);
50 } else {
51 if (a.shape() != shape) NDA_RUNTIME_ERROR << "Error in nda::detail::resize_or_check: Dimension mismatch: " << shape << " != " << a.shape();
52 }
53 }
54
55 // Given an array/view, prepare and return the corresponding h5::array_view to be written/read into.
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());
60 auto v = h5::array_interface::array_view{h5::hdf5_type<get_value_t<A>>(), (void *)a.data(), A::rank, is_complex_v<typename A::value_type>};
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];
65 }
66 return v;
67 }
68
69 // Create an h5::char_buf from an 1-dimensional nda::MemoryArray.
70 template <MemoryArrayOfRank<1> A>
71 h5::char_buf to_char_buf(A const &a) {
72 // get size of longest string
73 size_t s = 0;
74 for (auto &x : a) s = std::max(s, x.size() + 1);
75
76 // copy each string to a buffer and pad with zeros
77 std::vector<char> buf;
78 buf.resize(a.size() * s, 0x00);
79 size_t i = 0;
80 for (auto &x : a) {
81 strcpy(&buf[i * s], x.c_str());
82 ++i;
83 }
84
85 // return h5::char_buf
86 auto len = h5::v_t{size_t(a.size()), s};
87 return {buf, len};
88 }
89
90 // Write an h5::char_buf into a 1-dimensional nda::MemoryArray.
91 template <MemoryArrayOfRank<1> A>
92 void from_char_buf(h5::char_buf const &cb, A &a) {
93 // prepare the array
94 resize_or_check(a, std::array<long, 1>{static_cast<long>(cb.lengths[0])});
95
96 // loop over all strings
97 auto len_string = cb.lengths[1];
98 size_t i = 0;
99 for (auto &x : a) {
100 x = "";
101 x.append(&cb.buffer[i * len_string]);
102 ++i;
103 }
104 }
105
106 } // namespace detail
107
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>) {
126 // 1-dimensional array/view of strings
127 h5_write(g, name, detail::to_char_buf(a));
128 } else if constexpr (is_scalar_v<typename A::value_type>) {
129 // n-dimensional array/view of scalars
130 // make a copy if the array/view is not in C-order and write the copy
131 if (not a.indexmap().is_stride_order_C()) {
132 using h5_arr_t = nda::array<get_value_t<A>, A::rank>;
133 auto a_c_layout = h5_arr_t{a.shape()};
134 a_c_layout() = a;
135 h5_write(g, name, a_c_layout, compress);
136 return;
137 }
138
139 // prepare the h5::array_view and write it
140 auto v = detail::prepare_h5_array_view(a);
141 h5::array_interface::write(g, name, v, compress);
142 } else {
143 // n-dimensional array/view of some generic unknown type
144 auto g2 = g.create_group(name);
145 h5_write(g2, "shape", a.shape());
146 auto make_name = [](auto i0, auto... is) { return (std::to_string(i0) + ... + ("_" + std::to_string(is))); };
147 nda::for_each(a.shape(), [&](auto... is) { h5_write(g2, make_name(is...), a(is...)); });
148 }
149 }
150
166 template <size_t NDim, typename... IRs>
167 auto hyperslab_and_shape_from_slice(std::tuple<IRs...> const &slice, std::vector<h5::hsize_t> const &ds_shape, bool is_complex) {
168 // number of parameters specifying the slice
169 static constexpr auto size_of_slice = sizeof...(IRs);
170
171 // number of nda::ellipsis objects in the slice (only 1 is allowed)
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);
175
176 // position of the ellipsis in the slice parameters
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...>{});
181
182 // number of integers in the slice parameters
183 static constexpr auto integer_count = (std::integral<IRs> + ... + 0);
184
185 // number of nda::range and nda::range::all_t objects in the slice parameters
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");
189
190 // number of dimensions spanned by the ellipsis
191 static constexpr auto ellipsis_width = NDim - range_count;
192
193 // number of dimensions in the hyperslab
194 static constexpr auto slab_rank = NDim + integer_count;
195
196 // check rank of the dataset (must match the rank of the hyperslabs)
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;
201
202 // create and return the hyperslab and the shape array
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 {
206 (
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>) {
210 slab.offset[n] = ir;
211 slab.count[n] = 1;
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];
216 }
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();
222 } else {
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];
226 }
227 }(Is, std::get<Is>(slice)),
228 ...);
229 }(std::make_index_sequence<size_of_slice>{});
230 return std::make_pair(slab, shape);
231 }
232
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) {
251 // compile-time checks
252 constexpr int size_of_slice = sizeof...(IRs);
253 static_assert(size_of_slice > 0, "Error in nda::h5_write: Invalid empty slice");
254 static_assert(is_scalar_v<typename A::value_type>, "Error in nda::h5_write: Slicing is only supported for scalar types");
255 static constexpr bool is_complex = is_complex_v<typename A::value_type>;
256
257 // make a copy if the array/view is not in C-order and write the copy
258 if (not a.indexmap().is_stride_order_C()) {
259 using h5_arr_t = nda::array<get_value_t<A>, A::rank>;
260 auto a_c_layout = h5_arr_t{a.shape()};
261 a_c_layout() = a;
262 h5_write(g, name, a_c_layout, slice);
263 return;
264 }
265
266 // get dataset info and check that the datatypes of the dataset and the array match
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";
270
271 // get hyperslab and shape from the slice and check that the shapes match
272 auto const [sl, sh] = hyperslab_and_shape_from_slice<A::rank>(slice, ds_info.lengths, is_complex);
273 if (sh != a.shape()) NDA_RUNTIME_ERROR << "Error in nda::h5_write: Incompatible slice and array shape: " << sh << " != " << a.shape();
274
275 // prepare the h5::array_view and write it
276 auto v = detail::prepare_h5_array_view(a);
277 h5::array_interface::write_slice(g, name, v, sl);
278 }
279
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);
304
305 if constexpr (std::is_same_v<typename A::value_type, std::string>) {
306 // 1-dimensional array/view of strings
307 static_assert(!slicing, "Error in nda::h5_read: Slicing not supported for arrays/views of strings");
308 h5::char_buf cb;
309 h5_read(g, name, cb);
310 detail::from_char_buf(cb, a);
311 } else if constexpr (is_scalar_v<typename A::value_type>) {
312 // n-dimensional array/view of scalars
313 // read into a temporary array if the array/view is not in C-order and copy the elements
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());
319 a() = a_c_layout;
320 return;
321 }
322
323 // get dataset info
324 auto ds_info = h5::array_interface::get_dataset_info(g, name);
325
326 // allow to read non-complex data into a complex array
327 static constexpr bool is_complex = is_complex_v<typename A::value_type>;
328 if constexpr (is_complex) {
329 if (!ds_info.has_complex_attribute) {
331 h5_read(g, name, tmp, slice);
332 a = tmp;
333 return;
334 }
335 }
336
337 // get the hyperslab and the shape of the slice
338 std::array<long, A::rank> shape{};
339 auto slab = h5::array_interface::hyperslab{};
340 if constexpr (slicing) {
341 auto const [sl, sh] = hyperslab_and_shape_from_slice<A::rank>(slice, ds_info.lengths, is_complex);
342 slab = sl;
343 shape = sh;
344 } else {
345 for (int u = 0; u < A::rank; ++u) shape[u] = ds_info.lengths[u]; // NB : correct for complex
346 }
347
348 // resize the array or check that the shape matches
349 detail::resize_or_check(a, shape);
350
351 // check the rank of the dataset and the array
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;
355
356 // prepare the h5::array_view and read into it
357 auto v = detail::prepare_h5_array_view(a);
358 h5::array_interface::read(g, name, v, slab);
359 } else {
360 // n-dimensional array/view of some generic unknown type
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);
363
364 // get and check the shape or resize the array
365 std::array<long, A::rank> h5_shape;
366 h5_read(g2, "shape", h5_shape);
367 detail::resize_or_check(a, h5_shape);
368
369 // read element-by-element using the appropriate h5_read implementation
370 auto make_name = [](auto i0, auto... is) { return (std::to_string(i0) + ... + ("_" + std::to_string(is))); };
371 nda::for_each(a.shape(), [&](auto... is) { h5_read(g2, make_name(is...), a(is...)); });
372 }
373 }
374
376
377} // namespace nda
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.
Definition h5.hpp:302
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.
Definition h5.hpp:124
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,...
Definition h5.hpp:167
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.
Definition for_each.hpp:116
std::string to_string(std::array< T, R > const &a)
Get a string representation of a std::array.
Definition array.hpp:52
constexpr bool is_complex_v
Constexpr variable that is true if type T is a std::complex type.
Definition traits.hpp:65
constexpr bool is_scalar_v
Constexpr variable that is true if type S is a scalar type, i.e. arithmetic or complex.
Definition traits.hpp:69
Includes the itertools header and provides some additional utilities.
Provides type traits for the nda library.