Skip to content

Commit

Permalink
Add C++ docs and do some trivial clean up
Browse files Browse the repository at this point in the history
  • Loading branch information
@Thoemi09 @Wentzell
Thoemi09 authored and Wentzell committed Apr 24, 2024
1 parent f9da29e commit 0b446ec
Show file tree
Hide file tree
Showing 25 changed files with 1,797 additions and 888 deletions.
206 changes: 109 additions & 97 deletions c++/h5/array_interface.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,65 +15,88 @@
// Authors: Olivier Parcollet, Nils Wentzell

#include "./array_interface.hpp"
#include "./macros.hpp"
#include "./stl/string.hpp"

#include <hdf5.h>
#include <hdf5_hl.h>

#include <numeric>
#include <algorithm>
#include <iostream> // DEBUG
#include <iostream>
#include <stdexcept>
#include <utility>
#include <vector>

namespace h5::array_interface {

//------------------------------------------------
// the dataspace corresponding to the array. Contiguous data only...
dataspace make_mem_dspace(h5_array_view const &v) {
namespace {

if (v.rank() == 0) return H5Screate(H5S_SCALAR);
// Create an HDF5 memory dataspace.
dataspace make_mem_dspace(h5_array_view const &v) {
// scalar case
if (v.rank() == 0) return H5Screate(H5S_SCALAR);

dataspace dspace = H5Screate_simple(v.slab.rank(), v.L_tot.data(), nullptr);
if (!dspace.is_valid()) throw std::runtime_error("Cannot create the dataspace");
// create a dataspace of rank v.rank() and with shape v.L_tot
dataspace dspace = H5Screate_simple(v.slab.rank(), v.L_tot.data(), nullptr);
if (!dspace.is_valid()) throw std::runtime_error("Error in make_mem_dspace: Creating the dataspace for an h5_array_view failed");

herr_t err = H5Sselect_hyperslab(dspace, H5S_SELECT_SET, v.slab.offset.data(), v.slab.stride.data(), v.slab.count.data(),
(v.slab.block.empty() ? nullptr : v.slab.block.data()));
if (err < 0) throw std::runtime_error("Cannot set hyperslab");
// select the hyperslab according to v.slab
herr_t err = H5Sselect_hyperslab(dspace, H5S_SELECT_SET, v.slab.offset.data(), v.slab.stride.data(), v.slab.count.data(),
(v.slab.block.empty() ? nullptr : v.slab.block.data()));
if (err < 0) throw std::runtime_error("Error in make_mem_dspace: Selecting the hyperslab failed");

return dspace;
}
// return the dataspace
return dspace;
}

//------------------------------------------------
} // namespace

std::pair<v_t, v_t> get_L_tot_and_strides_h5(long const *stri, int rank, long total_size) {
std::pair<v_t, v_t> get_L_tot_and_strides_h5(long const *np_strides, int rank, long view_size) {
// scalar case: return empty vectors
if (rank == 0) return {};
if (total_size == 0) return {v_t(rank, 0), v_t(rank, 1)}; // if size = 0, we return (0,0,0), (1,1,1)

v_t Ltot(rank), strides_h5(rank);
for (int u = 0; u < rank; ++u) strides_h5[u] = stri[u];
Ltot[0] = total_size;
// empty nd-array case: return (0,0,0), (1,1,1)
if (view_size == 0) return {v_t(rank, 0), v_t(rank, 1)};

// general case
v_t Ltot(rank), h5_strides(rank);
for (int u = 0; u < rank; ++u) h5_strides[u] = np_strides[u];
Ltot[0] = view_size;

for (int u = rank - 2; u >= 0; --u) {
// L[u+1] as gcd of size and stride[u] ... stride[0]
long L = strides_h5[u];
// L becomes the gcd
for (int v = u - 1; v >= 0; --v) { L = std::gcd(L, strides_h5[v]); }
// L[u+1] as gcd of size and stride[u] ... stride[0]
hsize_t L = h5_strides[u];
// L becomes the gcd
for (int v = u - 1; v >= 0; --v) { L = std::gcd(L, h5_strides[v]); }
// divides
for (int v = u; v >= 0; --v) { strides_h5[v] /= L; }
for (int v = u; v >= 0; --v) { h5_strides[v] /= L; }
Ltot[u + 1] = L;
}

return {Ltot, strides_h5};
return {Ltot, h5_strides};
}

//-------------------------------------------------------
// write
//-------------------------------------------------------
h5_lengths_type get_h5_lengths_type(group g, std::string const &name) {
// open dataset
dataset ds = g.open_dataset(name);

void write(group g, std::string const &name, h5_array_view const &v, bool compress) {
// retrieve shape information
datatype ty = H5Dget_type(ds);
bool has_complex_attribute = H5LTfind_attribute(ds, "__complex__");
dataspace dspace = H5Dget_space(ds);
int rank = H5Sget_simple_extent_ndims(dspace);
v_t dims_out(rank);
H5Sget_simple_extent_dims(dspace, dims_out.data(), nullptr);

return {std::move(dims_out), ty, has_complex_attribute};
}

void write(group g, std::string const &name, h5_array_view const &v, bool compress) {
// unlink the dataset if it already exists
g.unlink(name);

// Some properties for the dataset : add compression
// chunk the dataset and add compression
proplist cparms = H5P_DEFAULT;
if (compress and (v.rank() != 0)) {
int n_dims = v.rank();
Expand All @@ -91,141 +114,130 @@ namespace h5::array_interface {
H5Pset_deflate(cparms, 1);
}

// dataspace for the dataset in the file
// dataspace for the dataset in the file: v.slab.block shape {1, ..., 1} is assumed
dataspace file_dspace = H5Screate_simple(v.slab.rank(), v.slab.count.data(), nullptr);

// create the dataset in the file
dataset ds = H5Dcreate2(g, name.c_str(), v.ty, file_dspace, H5P_DEFAULT, cparms, H5P_DEFAULT);
if (!ds.is_valid()) throw std::runtime_error("Cannot create the dataset " + name + " in the group " + g.name());
if (!ds.is_valid())
throw std::runtime_error("Error in h5::array_interface::write: Creating the dataset " + name + " in the group " + g.name() + " failed");

// memory data space
// memory dataspace
dataspace mem_dspace = make_mem_dspace(v);

// write to the file dataset
if (H5Sget_simple_extent_npoints(mem_dspace) > 0) { // avoid writing empty arrays
herr_t err = H5Dwrite(ds, v.ty, mem_dspace, H5S_ALL, H5P_DEFAULT, v.start);
if (err < 0) throw std::runtime_error("Error writing the scalar dataset " + name + " in the group" + g.name());
if (err < 0)
throw std::runtime_error("Error in h5::array_interface::write: Writing to the dataset " + name + " in the group" + g.name() + " failed");
}

// If we are dealing with complex, we had the attribute
// add complex attribute if the data is complex valued
if (v.is_complex) h5_write_attribute(ds, "__complex__", "1");
}

//-------------------------------------------------------

void write_slice(group g, std::string const &name, h5_array_view const &v, h5_lengths_type lt, hyperslab sl) {

// empty hyperslab
if (sl.empty()) return;

if (v.slab.count != sl.count) throw std::runtime_error("Error in h5::write_slice: Slab for memory and file incompatible");
// check consistency of input: block shape {1, ..., 1} is assumed
if (v.slab.count != sl.count) throw std::runtime_error("Error in h5::array_interface::write_slice: Memory and file slabs are incompatible");
if (not hdf5_type_equal(v.ty, lt.ty))
throw std::runtime_error("Error in h5::write_slice: Mismatching types. Expecting a " + get_name_of_h5_type(v.ty)
+ " while the array stored in the hdf5 file has type " + get_name_of_h5_type(lt.ty));
throw std::runtime_error("Error in h5::array_interface::write_slice: Incompatible HDF5 types: " + get_name_of_h5_type(v.ty)
+ " != " + get_name_of_h5_type(lt.ty));

// For a sliced write we assume the dataset already exists
// open existing dataset, get dataspace and select hyperslab
dataset ds = g.open_dataset(name);
dataspace file_dspace = H5Dget_space(ds);
herr_t err = H5Sselect_hyperslab(file_dspace, H5S_SELECT_SET, sl.offset.data(), sl.stride.data(), sl.count.data(),
(sl.block.empty() ? nullptr : sl.block.data()));
if (err < 0) throw std::runtime_error("Error in h5::array_interface::write_slice: Selecting the hyperslab failed");

herr_t err = H5Sselect_hyperslab(file_dspace, H5S_SELECT_SET, sl.offset.data(), sl.stride.data(), sl.count.data(),
(sl.block.empty() ? nullptr : sl.block.data()));
if (err < 0) throw std::runtime_error("Cannot set hyperslab");

// memory dataspace
dataspace mem_dspace = make_mem_dspace(v);

// write to the selected hyperslab of the file dataset
if (H5Sget_simple_extent_npoints(file_dspace) > 0) {
err = H5Dwrite(ds, v.ty, mem_dspace, file_dspace, H5P_DEFAULT, v.start);
if (err < 0) throw std::runtime_error("Error opening the scalar dataset " + name + " for sliced write in the group " + g.name());
if (err < 0)
throw std::runtime_error("Error in h5::array_interface::write_slice: Writing the dataset " + name + " in the group " + g.name() + " failed");
}
}

//-------------------------------------------------------------

void write_attribute(object obj, std::string const &name, h5_array_view v) {
// check if the attribute already exists
if (H5LTfind_attribute(obj, name.c_str()) != 0)
throw std::runtime_error("Error in h5::array_interface::write_attribute: Attribute " + name + " already exists");

if (H5LTfind_attribute(obj, name.c_str()) != 0) throw std::runtime_error("The attribute " + name + " is already present. Can not overwrite");

// memory dataspace
dataspace mem_dspace = make_mem_dspace(v);

// create attribute to write to
attribute attr = H5Acreate2(obj, name.c_str(), v.ty, mem_dspace, H5P_DEFAULT, H5P_DEFAULT);
if (!attr.is_valid()) throw std::runtime_error("Cannot create the attribute " + name);
if (!attr.is_valid()) throw std::runtime_error("Error in h5::array_interface::write_attribute: Creating the attribute " + name + " failed");

// write to the attribute
herr_t err = H5Awrite(attr, v.ty, v.start);
if (err < 0) throw std::runtime_error("Cannot write the attribute " + name);
if (err < 0) throw std::runtime_error("Error in h5::array_interface::write_attribute: Writing to the attribute " + name + " failed");
}

//-------------------------------------------------------
// READ
//-------------------------------------------------------

h5_lengths_type get_h5_lengths_type(group g, std::string const &name) {

dataset ds = g.open_dataset(name);

bool has_complex_attribute = H5LTfind_attribute(ds, "__complex__"); // the array in file should be interpreted as a complex
dataspace dspace = H5Dget_space(ds);
int rank = H5Sget_simple_extent_ndims(dspace);

// need to use hsize_t here and the vector is truncated at rank
v_t dims_out(rank);
H5Sget_simple_extent_dims(dspace, dims_out.data(), nullptr);

// get the type from the file
datatype ty = H5Dget_type(ds);
return {std::move(dims_out), ty, has_complex_attribute};
}

//--------------------------------------------------------

void read(group g, std::string const &name, h5_array_view v, h5_lengths_type lt, hyperslab sl) {

// open dataset and get dataspace
dataset ds = g.open_dataset(name);
dataspace file_dspace = H5Dget_space(ds);

// If provided, select the hyperslab of the file data
// if provided, select the hyperslab of the file dataspace
if (not sl.empty()) {
herr_t err = H5Sselect_hyperslab(file_dspace, H5S_SELECT_SET, sl.offset.data(), sl.stride.data(), sl.count.data(),
(sl.block.empty() ? nullptr : sl.block.data()));
if (err < 0) throw std::runtime_error("Cannot set hyperslab");
if (err < 0) throw std::runtime_error("Error in h5::array_interface::read: selecting the hyperslab failed");
}

// Checks
// check consistency of input
if (H5Tget_class(v.ty) != H5Tget_class(lt.ty))
throw std::runtime_error("Incompatible types in h5_read. Expecting a " + get_name_of_h5_type(v.ty)
+ " while the array stored in the hdf5 file has type " + get_name_of_h5_type(lt.ty));
throw std::runtime_error("Error in h5::array_interface::read: Incompatible HDF5 types: " + get_name_of_h5_type(v.ty)
+ " != " + get_name_of_h5_type(lt.ty));

if (not hdf5_type_equal(v.ty, lt.ty))
std::cerr << "WARNING: Mismatching types in h5_read. Expecting a " + get_name_of_h5_type(v.ty)
+ " while the array stored in the hdf5 file has type " + get_name_of_h5_type(lt.ty) + "\n";
std::cerr << "WARNING: HDF5 type mismatch while reading into an h5_array_view: " + get_name_of_h5_type(v.ty)
+ " != " + get_name_of_h5_type(lt.ty) + "\n";

if (lt.rank() != v.rank())
throw std::runtime_error("h5 read. Rank mismatch : expecting in file a rank " + std::to_string(v.rank())
+ " while the array stored in the hdf5 file has rank " + std::to_string(lt.rank()));
throw std::runtime_error("Error in h5::array_interface::read: Incompatible ranks: " + std::to_string(v.rank())
+ " != " + std::to_string(lt.rank()));

if (sl.empty() and lt.lengths != v.slab.count) throw std::runtime_error("h5 read. Lengths mismatch");
// block shape of {1, ..., 1} is assumed
if (sl.empty() and lt.lengths != v.slab.count) throw std::runtime_error("Error in h5::array_interface::read: Incompatible shapes");

// memory dataspace
dataspace mem_dspace = make_mem_dspace(v);

// read the selected hyperslab from the file dataset
if (H5Sget_simple_extent_npoints(file_dspace) > 0) {
herr_t err = H5Dread(ds, v.ty, mem_dspace, file_dspace, H5P_DEFAULT, v.start);
if (err < 0) throw std::runtime_error("Error reading the scalar dataset " + name + " in the group " + g.name());
if (err < 0)
throw std::runtime_error("Error in h5::array_interface::read: Reading the dataset " + name + " in the group " + g.name() + " failed");
}
}

//-------------------------------------------------------------

void read_attribute(object obj, std::string const &name, h5_array_view v) {

//if (v.rank() != 0) throw std::runtime_error("Non scalar attribute not implemented");

// open attribute
attribute attr = H5Aopen(obj, name.c_str(), H5P_DEFAULT);
if (!attr.is_valid()) throw std::runtime_error("Cannot open the attribute " + name);
if (!attr.is_valid()) throw std::runtime_error("Error in h5::array_interface::read_attribute: Opening the attribute " + name + " failed");

// get dataspace information
dataspace space = H5Aget_space(attr);
int rank = H5Sget_simple_extent_ndims(space);
if (rank != 0) throw std::runtime_error("Reading a scalar attribute and got rank !=0");
if (rank != 0) throw std::runtime_error("Error in h5::array_interface::read_attribute: Attribute " + name + " has a rank != 0");

// get datatype information
auto eq = H5Tequal(H5Aget_type(attr), v.ty);
if (eq < 0) throw std::runtime_error("Type comparison failure in reading attribute");
if (eq == 0) throw std::runtime_error("Type mismatch in reading attribute");
if (eq < 0) throw std::runtime_error("Error in h5::array_interface::read_attribute: H5Tequal call failed");
if (eq == 0) throw std::runtime_error("Error in h5::array_interface::read_attribute: Incompatible HDF5 types");

// read the attribute
auto err = H5Aread(attr, v.ty, v.start);
if (err < 0) throw std::runtime_error("Cannot read the attribute " + name);
if (err < 0) throw std::runtime_error("Error in h5::array_interface::read_attribute: Reading the attribute " + name + " failed");
}

} // namespace h5::array_interface
Loading

0 comments on commit 0b446ec

Please sign in to comment.