From 5ff6236a39929c9aa981554d8d061628a2795a84 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Wed, 18 Dec 2019 17:24:38 -0800 Subject: [PATCH] HDF5: Complex Support --- include/openPMD/IO/HDF5/HDF5Auxiliary.hpp | 181 +++++++++++------- include/openPMD/IO/HDF5/HDF5IOHandlerImpl.hpp | 4 + src/IO/HDF5/HDF5IOHandler.cpp | 123 ++++++++++-- 3 files changed, 220 insertions(+), 88 deletions(-) diff --git a/include/openPMD/IO/HDF5/HDF5Auxiliary.hpp b/include/openPMD/IO/HDF5/HDF5Auxiliary.hpp index 28445b47fb..a69063bf5f 100644 --- a/include/openPMD/IO/HDF5/HDF5Auxiliary.hpp +++ b/include/openPMD/IO/HDF5/HDF5Auxiliary.hpp @@ -27,82 +27,105 @@ #include +#include #include +#include +#include +#include +#include namespace openPMD { -inline hid_t -getH5DataType(Attribute const& att) +struct GetH5DataType { - using DT = Datatype; - switch( att.dtype ) + std::unordered_map< std::string, hid_t > m_userTypes; + + GetH5DataType( std::unordered_map< std::string, hid_t > userTypes ) + : m_userTypes{ std::move(userTypes) } { - case DT::CHAR: - case DT::VEC_CHAR: - return H5Tcopy(H5T_NATIVE_CHAR); - case DT::UCHAR: - case DT::VEC_UCHAR: - return H5Tcopy(H5T_NATIVE_UCHAR); - case DT::SHORT: - case DT::VEC_SHORT: - return H5Tcopy(H5T_NATIVE_SHORT); - case DT::INT: - case DT::VEC_INT: - return H5Tcopy(H5T_NATIVE_INT); - case DT::LONG: - case DT::VEC_LONG: - return H5Tcopy(H5T_NATIVE_LONG); - case DT::LONGLONG: - case DT::VEC_LONGLONG: - return H5Tcopy(H5T_NATIVE_LLONG); - case DT::USHORT: - case DT::VEC_USHORT: - return H5Tcopy(H5T_NATIVE_USHORT); - case DT::UINT: - case DT::VEC_UINT: - return H5Tcopy(H5T_NATIVE_UINT); - case DT::ULONG: - case DT::VEC_ULONG: - return H5Tcopy(H5T_NATIVE_ULONG); - case DT::ULONGLONG: - case DT::VEC_ULONGLONG: - return H5Tcopy(H5T_NATIVE_ULLONG); - case DT::FLOAT: - case DT::VEC_FLOAT: - return H5Tcopy(H5T_NATIVE_FLOAT); - case DT::DOUBLE: - case DT::ARR_DBL_7: - case DT::VEC_DOUBLE: - return H5Tcopy(H5T_NATIVE_DOUBLE); - case DT::LONG_DOUBLE: - case DT::VEC_LONG_DOUBLE: - return H5Tcopy(H5T_NATIVE_LDOUBLE); - case DT::STRING: - { - hid_t string_t_id = H5Tcopy(H5T_C_S1); - H5Tset_size(string_t_id, att.get< std::string >().size()); - return string_t_id; - } - case DT::VEC_STRING: - { - hid_t string_t_id = H5Tcopy(H5T_C_S1); - size_t max_len = 0; - for( std::string const& s : att.get< std::vector< std::string > >() ) - max_len = std::max(max_len, s.size()); - H5Tset_size(string_t_id, max_len); - return string_t_id; + } + + hid_t + operator()(Attribute const &att) + { + using DT = Datatype; + switch (att.dtype) { + case DT::CHAR: + case DT::VEC_CHAR: + return H5Tcopy(H5T_NATIVE_CHAR); + case DT::UCHAR: + case DT::VEC_UCHAR: + return H5Tcopy(H5T_NATIVE_UCHAR); + case DT::SHORT: + case DT::VEC_SHORT: + return H5Tcopy(H5T_NATIVE_SHORT); + case DT::INT: + case DT::VEC_INT: + return H5Tcopy(H5T_NATIVE_INT); + case DT::LONG: + case DT::VEC_LONG: + return H5Tcopy(H5T_NATIVE_LONG); + case DT::LONGLONG: + case DT::VEC_LONGLONG: + return H5Tcopy(H5T_NATIVE_LLONG); + case DT::USHORT: + case DT::VEC_USHORT: + return H5Tcopy(H5T_NATIVE_USHORT); + case DT::UINT: + case DT::VEC_UINT: + return H5Tcopy(H5T_NATIVE_UINT); + case DT::ULONG: + case DT::VEC_ULONG: + return H5Tcopy(H5T_NATIVE_ULONG); + case DT::ULONGLONG: + case DT::VEC_ULONGLONG: + return H5Tcopy(H5T_NATIVE_ULLONG); + case DT::FLOAT: + case DT::VEC_FLOAT: + return H5Tcopy(H5T_NATIVE_FLOAT); + case DT::DOUBLE: + case DT::ARR_DBL_7: + case DT::VEC_DOUBLE: + return H5Tcopy(H5T_NATIVE_DOUBLE); + case DT::LONG_DOUBLE: + case DT::VEC_LONG_DOUBLE: + return H5Tcopy(H5T_NATIVE_LDOUBLE); + + case DT::CFLOAT: + case DT::VEC_CFLOAT: + return H5Tcopy( m_userTypes.at( typeid(std::complex< float >).name() ) ); + case DT::CDOUBLE: + case DT::VEC_CDOUBLE: + return H5Tcopy( m_userTypes.at( typeid(std::complex< double >).name() ) ); + case DT::CLONG_DOUBLE: + case DT::VEC_CLONG_DOUBLE: + return H5Tcopy( m_userTypes.at( typeid(std::complex< long double >).name() ) ); + + case DT::STRING: { + hid_t string_t_id = H5Tcopy(H5T_C_S1); + H5Tset_size(string_t_id, att.get().size()); + return string_t_id; + } + case DT::VEC_STRING: { + hid_t string_t_id = H5Tcopy(H5T_C_S1); + size_t max_len = 0; + for (std::string const &s : att.get >()) + max_len = std::max(max_len, s.size()); + H5Tset_size(string_t_id, max_len); + return string_t_id; + } + case DT::BOOL: + return H5Tcopy( m_userTypes.at( typeid(bool).name() ) ); + case DT::DATATYPE: + throw std::runtime_error("[HDF5] Meta-Datatype leaked into IO"); + case DT::UNDEFINED: + throw std::runtime_error("[HDF5] Unknown Attribute datatype (HDF5 datatype)"); + default: + throw std::runtime_error("[HDF5] Datatype not implemented"); } - case DT::BOOL: - return H5Tcopy(H5T_NATIVE_HBOOL); - case DT::DATATYPE: - throw std::runtime_error("Meta-Datatype leaked into IO"); - case DT::UNDEFINED: - throw std::runtime_error("Unknown Attribute datatype (HDF5 datatype)"); - default: - throw std::runtime_error("Datatype not implemented in HDF5 IO"); } -} +}; inline hid_t getH5DataSpace(Attribute const& att) @@ -123,6 +146,9 @@ getH5DataSpace(Attribute const& att) case DT::FLOAT: case DT::DOUBLE: case DT::LONG_DOUBLE: + case DT::CFLOAT: + case DT::CDOUBLE: + case DT::CLONG_DOUBLE: case DT::STRING: case DT::BOOL: return H5Screate(H5S_SCALAR); @@ -217,6 +243,27 @@ getH5DataSpace(Attribute const& att) H5Sset_extent_simple(vec_t_id, 1, dims, nullptr); return vec_t_id; } + case DT::VEC_CFLOAT: + { + hid_t vec_t_id = H5Screate(H5S_SIMPLE); + hsize_t dims[1] = {att.get< std::vector< std::complex< float > > >().size()}; + H5Sset_extent_simple(vec_t_id, 1, dims, nullptr); + return vec_t_id; + } + case DT::VEC_CDOUBLE: + { + hid_t vec_t_id = H5Screate(H5S_SIMPLE); + hsize_t dims[1] = {att.get< std::vector< std::complex< double > > >().size()}; + H5Sset_extent_simple(vec_t_id, 1, dims, nullptr); + return vec_t_id; + } + case DT::VEC_CLONG_DOUBLE: + { + hid_t vec_t_id = H5Screate(H5S_SIMPLE); + hsize_t dims[1] = {att.get< std::vector< std::complex< long double > > >().size()}; + H5Sset_extent_simple(vec_t_id, 1, dims, nullptr); + return vec_t_id; + } case DT::VEC_STRING: { hid_t vec_t_id = H5Screate(H5S_SIMPLE); diff --git a/include/openPMD/IO/HDF5/HDF5IOHandlerImpl.hpp b/include/openPMD/IO/HDF5/HDF5IOHandlerImpl.hpp index ee0e2b0c10..f405d10e53 100644 --- a/include/openPMD/IO/HDF5/HDF5IOHandlerImpl.hpp +++ b/include/openPMD/IO/HDF5/HDF5IOHandlerImpl.hpp @@ -63,7 +63,11 @@ namespace openPMD hid_t m_datasetTransferProperty; hid_t m_fileAccessProperty; + // h5py compatible types for bool and complex hid_t m_H5T_BOOL_ENUM; + hid_t m_H5T_CFLOAT; + hid_t m_H5T_CDOUBLE; + hid_t m_H5T_CLONG_DOUBLE; }; // HDF5IOHandlerImpl #else class HDF5IOHandlerImpl diff --git a/src/IO/HDF5/HDF5IOHandler.cpp b/src/IO/HDF5/HDF5IOHandler.cpp index 6dda17b41a..091e35c6bb 100644 --- a/src/IO/HDF5/HDF5IOHandler.cpp +++ b/src/IO/HDF5/HDF5IOHandler.cpp @@ -30,10 +30,14 @@ # include "openPMD/IO/HDF5/HDF5FilePosition.hpp" #endif +#include #include #include #include #include +#include +#include +#include #include namespace openPMD @@ -49,9 +53,13 @@ HDF5IOHandlerImpl::HDF5IOHandlerImpl(AbstractIOHandler* handler) : AbstractIOHandlerImpl(handler), m_datasetTransferProperty{H5P_DEFAULT}, m_fileAccessProperty{H5P_DEFAULT}, - m_H5T_BOOL_ENUM{H5Tenum_create(H5T_NATIVE_INT8)} + m_H5T_BOOL_ENUM{H5Tenum_create(H5T_NATIVE_INT8)}, + m_H5T_CFLOAT{H5Tcreate(H5T_COMPOUND, sizeof(float) * 2)}, + m_H5T_CDOUBLE{H5Tcreate(H5T_COMPOUND, sizeof(double) * 2)}, + m_H5T_CLONG_DOUBLE{H5Tcreate(H5T_COMPOUND, sizeof(long double) * 2)} { - VERIFY(m_H5T_BOOL_ENUM >= 0, "[HDF5] Internal error: Failed to create HDF5 enum"); + // create a h5py compatible bool type + VERIFY(m_H5T_BOOL_ENUM >= 0, "[HDF5] Internal error: Failed to create bool enum"); std::string t{"TRUE"}; std::string f{"FALSE"}; int64_t tVal = 1; @@ -61,6 +69,17 @@ HDF5IOHandlerImpl::HDF5IOHandlerImpl(AbstractIOHandler* handler) VERIFY(status == 0, "[HDF5] Internal error: Failed to insert into HDF5 enum"); status = H5Tenum_insert(m_H5T_BOOL_ENUM, f.c_str(), &fVal); VERIFY(status == 0, "[HDF5] Internal error: Failed to insert into HDF5 enum"); + + // create h5py compatible complex types + VERIFY(m_H5T_CFLOAT >= 0, "[HDF5] Internal error: Failed to create complex float"); + VERIFY(m_H5T_CDOUBLE >= 0, "[HDF5] Internal error: Failed to create complex double"); + VERIFY(m_H5T_CLONG_DOUBLE >= 0, "[HDF5] Internal error: Failed to create complex long double"); + H5Tinsert(m_H5T_CFLOAT, "r", 0, H5T_NATIVE_FLOAT); + H5Tinsert(m_H5T_CFLOAT, "i", sizeof(float), H5T_NATIVE_FLOAT); + H5Tinsert(m_H5T_CDOUBLE, "r", 0, H5T_NATIVE_DOUBLE); + H5Tinsert(m_H5T_CDOUBLE, "i", sizeof(double), H5T_NATIVE_DOUBLE); + H5Tinsert(m_H5T_CLONG_DOUBLE, "r", 0, H5T_NATIVE_LDOUBLE); + H5Tinsert(m_H5T_CLONG_DOUBLE, "i", sizeof(long double), H5T_NATIVE_LDOUBLE); } HDF5IOHandlerImpl::~HDF5IOHandlerImpl() @@ -68,26 +87,36 @@ HDF5IOHandlerImpl::~HDF5IOHandlerImpl() herr_t status; status = H5Tclose(m_H5T_BOOL_ENUM); if( status < 0 ) - std::cerr << "Internal error: Failed to close HDF5 enum\n"; + std::cerr << "[HDF5] Internal error: Failed to close bool enum\n"; + status = H5Tclose(m_H5T_CFLOAT); + if( status < 0 ) + std::cerr << "[HDF5] Internal error: Failed to close complex float type\n"; + status = H5Tclose(m_H5T_CDOUBLE); + if( status < 0 ) + std::cerr << "[HDF5] Internal error: Failed to close complex double type\n"; + status = H5Tclose(m_H5T_CLONG_DOUBLE); + if( status < 0 ) + std::cerr << "[HDF5] Internal error: Failed to close complex long double type\n"; + while( !m_openFileIDs.empty() ) { auto file = m_openFileIDs.begin(); status = H5Fclose(*file); if( status < 0 ) - std::cerr << "Internal error: Failed to close HDF5 file (serial)\n"; + std::cerr << "[HDF5] Internal error: Failed to close HDF5 file (serial)\n"; m_openFileIDs.erase(file); } if( m_datasetTransferProperty != H5P_DEFAULT ) { status = H5Pclose(m_datasetTransferProperty); if( status < 0 ) - std::cerr << "Internal error: Failed to close HDF5 dataset transfer property\n"; + std::cerr << "[HDF5] Internal error: Failed to close HDF5 dataset transfer property\n"; } if( m_fileAccessProperty != H5P_DEFAULT ) { status = H5Pclose(m_fileAccessProperty); if( status < 0 ) - std::cerr << "Internal error: Failed to close HDF5 file access property\n"; + std::cerr << "[HDF5] Internal error: Failed to close HDF5 file access property\n"; } } @@ -215,7 +244,7 @@ HDF5IOHandlerImpl::createDataset(Writable* writable, if( d == Datatype::UNDEFINED ) { // TODO handle unknown dtype - std::cerr << "Datatype::UNDEFINED caught during dataset creation (serial HDF5)" << std::endl; + std::cerr << "[HDF5] Datatype::UNDEFINED caught during dataset creation (serial HDF5)" << std::endl; d = Datatype::BOOL; } Attribute a(0); @@ -239,7 +268,7 @@ HDF5IOHandlerImpl::createDataset(Writable* writable, std::string const& compression = parameters.compression; if( !compression.empty() ) - std::cerr << "Compression not yet implemented in HDF5 backend." + std::cerr << "[HDF5] Compression not yet implemented in HDF5 backend." << std::endl; /* { @@ -251,11 +280,11 @@ HDF5IOHandlerImpl::createDataset(Writable* writable, status = H5Pset_deflate(datasetCreationProperty, std::stoi(args[1])); VERIFY(status == 0, "[HDF5] Internal error: Failed to set deflate compression during dataset creation"); } else if( format == "szip" || format == "nbit" || format == "scaleoffset" ) - std::cerr << "Compression format " << format + std::cerr << "[HDF5] Compression format " << format << " not yet implemented. Data will not be compressed!" << std::endl; else - std::cerr << "Compression format " << format + std::cerr << "[HDF5] Compression format " << format << " unknown. Data will not be compressed!" << std::endl; } @@ -263,9 +292,15 @@ HDF5IOHandlerImpl::createDataset(Writable* writable, std::string const& transform = parameters.transform; if( !transform.empty() ) - std::cerr << "Custom transform not yet implemented in HDF5 backend." + std::cerr << "[HDF5] Custom transform not yet implemented in HDF5 backend." << std::endl; + GetH5DataType getH5DataType({ + { typeid(bool).name(), m_H5T_BOOL_ENUM }, + { typeid(std::complex< float >).name(), m_H5T_CFLOAT }, + { typeid(std::complex< double >).name(), m_H5T_CDOUBLE }, + { typeid(std::complex< long double >).name(), m_H5T_CLONG_DOUBLE }, + }); hid_t datatype = getH5DataType(a); VERIFY(datatype >= 0, "[HDF5] Internal error: Failed to get HDF5 datatype during dataset creation"); hid_t group_id = H5Dcreate(node_id, @@ -460,6 +495,14 @@ HDF5IOHandlerImpl::openDataset(Writable* writable, d = DT::FLOAT; else if( H5Tequal(dataset_type, H5T_NATIVE_DOUBLE) ) d = DT::DOUBLE; + else if( H5Tequal(dataset_type, H5T_NATIVE_LDOUBLE) ) + d = DT::LONG_DOUBLE; + else if( H5Tequal(dataset_type, m_H5T_CFLOAT) ) + d = DT::CFLOAT; + else if( H5Tequal(dataset_type, m_H5T_CDOUBLE) ) + d = DT::CDOUBLE; + else if( H5Tequal(dataset_type, m_H5T_CLONG_DOUBLE) ) + d = DT::CLONG_DOUBLE; else if( H5Tequal(dataset_type, H5T_NATIVE_USHORT) ) d = DT::USHORT; else if( H5Tequal(dataset_type, H5T_NATIVE_UINT) ) @@ -692,6 +735,13 @@ HDF5IOHandlerImpl::writeDataset(Writable* writable, std::shared_ptr< void const > data = parameters.data; + GetH5DataType getH5DataType({ + { typeid(bool).name(), m_H5T_BOOL_ENUM }, + { typeid(std::complex< float >).name(), m_H5T_CFLOAT }, + { typeid(std::complex< double >).name(), m_H5T_CDOUBLE }, + { typeid(std::complex< long double >).name(), m_H5T_CLONG_DOUBLE }, + }); + //TODO Check if parameter dtype and dataset dtype match Attribute a(0); a.dtype = parameters.dtype; @@ -700,8 +750,12 @@ HDF5IOHandlerImpl::writeDataset(Writable* writable, switch( a.dtype ) { using DT = Datatype; + case DT::LONG_DOUBLE: case DT::DOUBLE: case DT::FLOAT: + case DT::CLONG_DOUBLE: + case DT::CDOUBLE: + case DT::CFLOAT: case DT::SHORT: case DT::INT: case DT::LONG: @@ -758,11 +812,13 @@ HDF5IOHandlerImpl::writeAttribute(Writable* writable, Attribute const att(parameters.resource); Datatype dtype = parameters.dtype; herr_t status; - hid_t dataType; - if( dtype == Datatype::BOOL ) - dataType = m_H5T_BOOL_ENUM; - else - dataType = getH5DataType(att); + GetH5DataType getH5DataType({ + { typeid(bool).name(), m_H5T_BOOL_ENUM }, + { typeid(std::complex< float >).name(), m_H5T_CFLOAT }, + { typeid(std::complex< double >).name(), m_H5T_CDOUBLE }, + { typeid(std::complex< long double >).name(), m_H5T_CLONG_DOUBLE }, + }); + hid_t dataType = getH5DataType(att); VERIFY(dataType >= 0, "[HDF5] Internal error: Failed to get HDF5 datatype during attribute write"); std::string name = parameters.name; if( H5Aexists(node_id, name.c_str()) == 0 ) @@ -867,6 +923,24 @@ HDF5IOHandlerImpl::writeAttribute(Writable* writable, status = H5Awrite(attribute_id, dataType, &d); break; } + case DT::CFLOAT: + { + std::complex< float > f = att.get< std::complex< float > >(); + status = H5Awrite(attribute_id, dataType, &f); + break; + } + case DT::CDOUBLE: + { + std::complex< double > d = att.get< std::complex< double > >(); + status = H5Awrite(attribute_id, dataType, &d); + break; + } + case DT::CLONG_DOUBLE: + { + std::complex< long double > d = att.get< std::complex< long double > >(); + status = H5Awrite(attribute_id, dataType, &d); + break; + } case DT::STRING: status = H5Awrite(attribute_id, dataType, @@ -968,11 +1042,8 @@ HDF5IOHandlerImpl::writeAttribute(Writable* writable, } VERIFY(status == 0, "[HDF5] Internal error: Failed to write attribute " + name + " at " + concrete_h5_file_position(writable)); - if( dataType != m_H5T_BOOL_ENUM ) - { - status = H5Tclose(dataType); - VERIFY(status == 0, "[HDF5] Internal error: Failed to close HDF5 datatype during Attribute write"); - } + status = H5Tclose(dataType); + VERIFY(status == 0, "[HDF5] Internal error: Failed to close HDF5 datatype during Attribute write"); status = H5Aclose(attribute_id); VERIFY(status == 0, "[HDF5] Internal error: Failed to close attribute " + name + " at " + concrete_h5_file_position(writable) + " during attribute write"); @@ -1021,8 +1092,12 @@ HDF5IOHandlerImpl::readDataset(Writable* writable, switch( a.dtype ) { using DT = Datatype; + case DT::LONG_DOUBLE: case DT::DOUBLE: case DT::FLOAT: + case DT::CLONG_DOUBLE: + case DT::CDOUBLE: + case DT::CFLOAT: case DT::SHORT: case DT::INT: case DT::LONG: @@ -1042,6 +1117,12 @@ HDF5IOHandlerImpl::readDataset(Writable* writable, default: throw std::runtime_error("[HDF5] Datatype not implemented in HDF5 IO"); } + GetH5DataType getH5DataType({ + { typeid(bool).name(), m_H5T_BOOL_ENUM }, + { typeid(std::complex< float >).name(), m_H5T_CFLOAT }, + { typeid(std::complex< double >).name(), m_H5T_CDOUBLE }, + { typeid(std::complex< long double >).name(), m_H5T_CLONG_DOUBLE }, + }); hid_t dataType = getH5DataType(a); VERIFY(dataType >= 0, "[HDF5] Internal error: Failed to get HDF5 datatype during dataset read"); status = H5Dread(dataset_id,