From fe3c1737e84334cff6172b5db9df22077c22a180 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Thu, 3 Sep 2020 07:45:51 -0700 Subject: [PATCH] Support for Complex Types (#639) Add support for complex numbers as fundamental data type. Close #203 - [x] Tests - [x] internals: new types, comparisons, allowed casts, record components, patches and attributes, etc. - [x] JSON (writes to pairs, read re-identification as complex is todo) - [x] HDF5 - [x] ADIOS1 (no long double complex; no attributes of arrays of complex: https://github.com/ornladios/ADIOS/issues/212) - [x] ADIOS2: (complex attributes: https://github.com/ornladios/ADIOS2/issues/1908); (no long double complex https://github.com/ornladios/ADIOS2/pull/1907) - [x] Python bindings --- CHANGELOG.rst | 1 + docs/source/backends/adios1.rst | 4 + include/openPMD/Datatype.hpp | 138 ++++++++++++- include/openPMD/IO/ADIOS/ADIOS1Auxiliary.hpp | 9 + include/openPMD/IO/ADIOS/ADIOS2Auxiliary.hpp | 22 ++ include/openPMD/IO/ADIOS/ADIOS2IOHandler.hpp | 53 ++++- include/openPMD/IO/HDF5/HDF5Auxiliary.hpp | 179 ++++++++++------ include/openPMD/IO/HDF5/HDF5IOHandlerImpl.hpp | 4 + include/openPMD/IO/JSON/JSONIOHandlerImpl.hpp | 11 + include/openPMD/RecordComponent.hpp | 4 +- include/openPMD/auxiliary/Memory.hpp | 16 ++ include/openPMD/backend/Attribute.hpp | 17 ++ include/openPMD/binding/python/Numpy.hpp | 30 ++- src/Datatype.cpp | 50 ++++- src/IO/ADIOS/ADIOS2Auxiliary.cpp | 3 + src/IO/ADIOS/ADIOS2IOHandler.cpp | 5 +- src/IO/ADIOS/CommonADIOS1IOHandler.cpp | 61 +++++- src/IO/HDF5/HDF5IOHandler.cpp | 193 ++++++++++++++++-- src/IO/JSON/JSONIOHandlerImpl.cpp | 3 + src/RecordComponent.cpp | 10 + src/backend/Attributable.cpp | 19 ++ src/binding/python/Attributable.cpp | 33 ++- src/binding/python/RecordComponent.cpp | 16 ++ test/CoreTest.cpp | 7 + test/SerialIOTest.cpp | 77 +++++++ test/python/unittest/API/APITest.py | 92 ++++++++- 26 files changed, 948 insertions(+), 109 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 1e68d1bce4..40e27ad4e4 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -19,6 +19,7 @@ Features - ``Record(Component)``: ``scalar()``, ``constant()``, ``empty()`` #711 - Advanced backend configuration via JSON #569 #733 +- Support for complex floating point types #639 - Functionality to close an iteration (and associated files) #746 - Python: diff --git a/docs/source/backends/adios1.rst b/docs/source/backends/adios1.rst index bb2085143a..6bdebfcd52 100644 --- a/docs/source/backends/adios1.rst +++ b/docs/source/backends/adios1.rst @@ -73,6 +73,10 @@ Limitations Mea culpa, we did better in the past (in PIConGPU). Please consider using our ADIOS2 backend instead, on which we focus our developments these days. +.. note:: + + ADIOS1 does not support attributes that are `arrays of complex types `_. + Selected References ------------------- diff --git a/include/openPMD/Datatype.hpp b/include/openPMD/Datatype.hpp index 308d20fab1..6dd3b15085 100644 --- a/include/openPMD/Datatype.hpp +++ b/include/openPMD/Datatype.hpp @@ -1,4 +1,4 @@ -/* Copyright 2017-2020 Fabian Koller and Franz Poeschel +/* Copyright 2017-2020 Fabian Koller, Franz Poeschel, Axel Huebl * * This file is part of openPMD-api. * @@ -31,6 +31,8 @@ #include #include #include +#include + namespace openPMD @@ -43,6 +45,7 @@ enum class Datatype : int SHORT, INT, LONG, LONGLONG, USHORT, UINT, ULONG, ULONGLONG, FLOAT, DOUBLE, LONG_DOUBLE, + CFLOAT, CDOUBLE, CLONG_DOUBLE, STRING, VEC_CHAR, VEC_SHORT, @@ -57,6 +60,9 @@ enum class Datatype : int VEC_FLOAT, VEC_DOUBLE, VEC_LONG_DOUBLE, + VEC_CFLOAT, + VEC_CDOUBLE, + VEC_CLONG_DOUBLE, VEC_STRING, ARR_DBL_7, @@ -136,6 +142,9 @@ determineDatatype() else if( decay_equiv< T, float >::value ){ return DT::FLOAT; } else if( decay_equiv< T, double >::value ){ return DT::DOUBLE; } else if( decay_equiv< T, long double >::value ){ return DT::LONG_DOUBLE; } + else if( decay_equiv< T, std::complex< float > >::value ){ return DT::CFLOAT; } + else if( decay_equiv< T, std::complex< double > >::value ){ return DT::CDOUBLE; } + else if( decay_equiv< T, std::complex< long double > >::value ){ return DT::CLONG_DOUBLE; } else if( decay_equiv< T, std::string >::value ){ return DT::STRING; } else if( decay_equiv< T, std::vector< char > >::value ){ return DT::VEC_CHAR; } else if( decay_equiv< T, std::vector< short > >::value ){ return DT::VEC_SHORT; } @@ -150,6 +159,9 @@ determineDatatype() else if( decay_equiv< T, std::vector< float > >::value ){ return DT::VEC_FLOAT; } else if( decay_equiv< T, std::vector< double > >::value ){ return DT::VEC_DOUBLE; } else if( decay_equiv< T, std::vector< long double > >::value ){ return DT::VEC_LONG_DOUBLE; } + else if( decay_equiv< T, std::vector< std::complex< float > > >::value ){ return DT::VEC_CFLOAT; } + else if( decay_equiv< T, std::vector< std::complex< double > > >::value ){ return DT::VEC_CDOUBLE; } + else if( decay_equiv< T, std::vector< std::complex< long double > > >::value ){ return DT::VEC_CLONG_DOUBLE; } else if( decay_equiv< T, std::vector< std::string > >::value ){ return DT::VEC_STRING; } else if( decay_equiv< T, std::array< double, 7 > >::value ){ return DT::ARR_DBL_7; } else if( decay_equiv< T, bool >::value ){ return DT::BOOL; } @@ -178,6 +190,9 @@ determineDatatype(std::shared_ptr< T >) else if( decay_equiv< T, float >::value ){ return DT::FLOAT; } else if( decay_equiv< T, double >::value ){ return DT::DOUBLE; } else if( decay_equiv< T, long double >::value ){ return DT::LONG_DOUBLE; } + else if( decay_equiv< T, std::complex< float > >::value ){ return DT::CFLOAT; } + else if( decay_equiv< T, std::complex< double > >::value ){ return DT::CDOUBLE; } + else if( decay_equiv< T, std::complex< long double > >::value ){ return DT::CLONG_DOUBLE; } else if( decay_equiv< T, std::string >::value ){ return DT::STRING; } else if( decay_equiv< T, std::vector< char > >::value ){ return DT::VEC_CHAR; } else if( decay_equiv< T, std::vector< short > >::value ){ return DT::VEC_SHORT; } @@ -192,6 +207,9 @@ determineDatatype(std::shared_ptr< T >) else if( decay_equiv< T, std::vector< float > >::value ){ return DT::VEC_FLOAT; } else if( decay_equiv< T, std::vector< double > >::value ){ return DT::VEC_DOUBLE; } else if( decay_equiv< T, std::vector< long double > >::value ){ return DT::VEC_LONG_DOUBLE; } + else if( decay_equiv< T, std::vector< std::complex< float > > >::value ){ return DT::VEC_CFLOAT; } + else if( decay_equiv< T, std::vector< std::complex< double > > >::value ){ return DT::VEC_CDOUBLE; } + else if( decay_equiv< T, std::vector< std::complex< long double > > >::value ){ return DT::VEC_CLONG_DOUBLE; } else if( decay_equiv< T, std::vector< std::string > >::value ){ return DT::VEC_STRING; } else if( decay_equiv< T, std::array< double, 7 > >::value ){ return DT::ARR_DBL_7; } else if( decay_equiv< T, bool >::value ){ return DT::BOOL; } @@ -255,6 +273,19 @@ toBytes( Datatype d ) case DT::LONG_DOUBLE: case DT::VEC_LONG_DOUBLE: return sizeof(long double); + break; + case DT::CFLOAT: + case DT::VEC_CFLOAT: + return sizeof(float) * 2; + break; + case DT::CDOUBLE: + case DT::VEC_CDOUBLE: + return sizeof(double) * 2; + break; + case DT::CLONG_DOUBLE: + case DT::VEC_CLONG_DOUBLE: + return sizeof(long double) * 2; + break; case DT::BOOL: return sizeof(bool); case DT::DATATYPE: @@ -300,6 +331,9 @@ isVector( Datatype d ) case DT::VEC_FLOAT: case DT::VEC_DOUBLE: case DT::VEC_LONG_DOUBLE: + case DT::VEC_CFLOAT: + case DT::VEC_CDOUBLE: + case DT::VEC_CLONG_DOUBLE: case DT::VEC_STRING: return true; default: @@ -327,6 +361,35 @@ isFloatingPoint( Datatype d ) case DT::VEC_DOUBLE: case DT::LONG_DOUBLE: case DT::VEC_LONG_DOUBLE: + // note: complex floats are not std::is_floating_point + return true; + break; + default: + return false; + break; + } +} + +/** Compare if a Datatype is a complex floating point type + * + * Includes our vector types + * + * @param d Datatype to test + * @return true if complex floating point, otherwise false + */ +inline bool +isComplexFloatingPoint( Datatype d ) +{ + using DT = Datatype; + + switch( d ) + { + case DT::CFLOAT: + case DT::VEC_CFLOAT: + case DT::CDOUBLE: + case DT::VEC_CDOUBLE: + case DT::CLONG_DOUBLE: + case DT::VEC_CLONG_DOUBLE: return true; default: return false; @@ -349,6 +412,22 @@ isFloatingPoint() return isFloatingPoint( dtype ); } +/** Compare if a type is a complex floating point type + * + * Like isFloatingPoint but for complex floats + * + * @tparam T type to test + * @return true if complex floating point, otherwise false + */ +template< typename T > +inline bool +isComplexFloatingPoint() +{ + Datatype dtype = determineDatatype< T >(); + + return isComplexFloatingPoint(dtype); +} + /** Compare if a Datatype is an integer type * * contrary to std::is_integer, the types bool and char types are not @@ -430,6 +509,32 @@ isSameFloatingPoint( Datatype d ) return false; } +/** Compare if a Datatype is equivalent to a complex floating point type + * + * @tparam T_CFP complex floating point type to compare + * @param d Datatype to compare + * @return true if both types are complex floating point and same bitness, else false + */ +template< typename T_CFP > +inline bool +isSameComplexFloatingPoint( Datatype d ) +{ + // template + bool tt_is_cfp = isComplexFloatingPoint< T_CFP >(); + + // Datatype + bool dt_is_cfp = isComplexFloatingPoint( d ); + + if( + tt_is_cfp && + dt_is_cfp && + toBits( d ) == toBits( determineDatatype< T_CFP >() ) + ) + return true; + else + return false; +} + /** Compare if a Datatype is equivalent to an integer type * * @tparam T_Int signed or unsigned integer type to compare @@ -501,6 +606,18 @@ isSame( openPMD::Datatype const d, openPMD::Datatype const e ) ) return true; + // same complex floating point + bool d_is_cfp = isComplexFloatingPoint(d); + bool e_is_cfp = isComplexFloatingPoint(e); + + if( + d_is_cfp && + e_is_cfp && + d_is_vec == e_is_vec && + toBits( d ) == toBits( e ) + ) + return true; + return false; } @@ -571,6 +688,15 @@ ReturnType switchType( Datatype dt, Action action, Args &&... args ) case Datatype::LONG_DOUBLE: return action.OPENPMD_TEMPLATE_OPERATOR( )< long double >( std::forward< Args >( args )... ); + case Datatype::CFLOAT: + return action.OPENPMD_TEMPLATE_OPERATOR( )< std::complex< float > >( + std::forward< Args >( args )... ); + case Datatype::CDOUBLE: + return action.OPENPMD_TEMPLATE_OPERATOR( )< std::complex< double > >( + std::forward< Args >( args )... ); + case Datatype::CLONG_DOUBLE: + return action.OPENPMD_TEMPLATE_OPERATOR( )< std::complex< long double > >( + std::forward< Args >( args )... ); case Datatype::STRING: return action.OPENPMD_TEMPLATE_OPERATOR( )< std::string >( std::forward< Args >( args )... ); @@ -619,6 +745,16 @@ ReturnType switchType( Datatype dt, Action action, Args &&... args ) return action .OPENPMD_TEMPLATE_OPERATOR( )< std::vector< long double > >( std::forward< Args >( args )... ); + case Datatype::VEC_CFLOAT: + return action.OPENPMD_TEMPLATE_OPERATOR( )< std::vector< std::complex< float > > >( + std::forward< Args >( args )... ); + case Datatype::VEC_CDOUBLE: + return action.OPENPMD_TEMPLATE_OPERATOR( )< std::vector< std::complex< double > > >( + std::forward< Args >( args )... ); + case Datatype::VEC_CLONG_DOUBLE: + return action + .OPENPMD_TEMPLATE_OPERATOR( )< std::vector< std::complex< long double > > >( + std::forward< Args >( args )... ); case Datatype::VEC_STRING: return action .OPENPMD_TEMPLATE_OPERATOR( )< std::vector< std::string > >( diff --git a/include/openPMD/IO/ADIOS/ADIOS1Auxiliary.hpp b/include/openPMD/IO/ADIOS/ADIOS1Auxiliary.hpp index 389fd66b27..d18f2e5535 100644 --- a/include/openPMD/IO/ADIOS/ADIOS1Auxiliary.hpp +++ b/include/openPMD/IO/ADIOS/ADIOS1Auxiliary.hpp @@ -181,6 +181,15 @@ getBP1DataType(Datatype dtype) case DT::LONG_DOUBLE: case DT::VEC_LONG_DOUBLE: return adios_long_double; + case DT::CFLOAT: + case DT::VEC_CFLOAT: + return adios_complex; + case DT::CDOUBLE: + case DT::VEC_CDOUBLE: + return adios_double_complex; + case DT::CLONG_DOUBLE: + case DT::VEC_CLONG_DOUBLE: + throw unsupported_data_error("No native equivalent for Datatype::CLONG_DOUBLE found."); case DT::STRING: return adios_string; case DT::VEC_STRING: diff --git a/include/openPMD/IO/ADIOS/ADIOS2Auxiliary.hpp b/include/openPMD/IO/ADIOS/ADIOS2Auxiliary.hpp index 57411a21d3..a3c7b6d13a 100644 --- a/include/openPMD/IO/ADIOS/ADIOS2Auxiliary.hpp +++ b/include/openPMD/IO/ADIOS/ADIOS2Auxiliary.hpp @@ -25,6 +25,8 @@ #if openPMD_HAVE_ADIOS2 #include "openPMD/Datatype.hpp" #include +#include +#include #include #include @@ -79,12 +81,32 @@ namespace detail getSize( adios2::IO &, std::string const & attributeName ); }; + template < > struct AttributeInfoHelper< std::complex< long double > > + { + static typename std::vector< long double >::size_type + getSize( adios2::IO &, std::string const & ) + { + throw std::runtime_error( + "[ADIOS2] Internal error: no support for long double complex attribute types" ); + } + }; + template < typename T > struct AttributeInfoHelper< std::vector< T > > { static typename std::vector< T >::size_type getSize( adios2::IO &, std::string const & attributeName ); }; + template < > struct AttributeInfoHelper< std::vector< std::complex< long double > > > + { + static typename std::vector< std::complex< long double > >::size_type + getSize( adios2::IO &, std::string const & ) + { + throw std::runtime_error( + "[ADIOS2] Internal error: no support for long double complex vector attribute types" ); + } + }; + template < typename T, std::size_t n > struct AttributeInfoHelper< std::array< T, n > > { diff --git a/include/openPMD/IO/ADIOS/ADIOS2IOHandler.hpp b/include/openPMD/IO/ADIOS/ADIOS2IOHandler.hpp index ccaf7afe26..0ef7085f50 100644 --- a/include/openPMD/IO/ADIOS/ADIOS2IOHandler.hpp +++ b/include/openPMD/IO/ADIOS/ADIOS2IOHandler.hpp @@ -430,6 +430,48 @@ namespace detail std::shared_ptr< Attribute::resource > resource ); }; + template< > struct AttributeTypes< std::complex< long double > > + { + using Attr = adios2::Attribute< std::complex< double > >; + using BasicType = double; + + static Attr createAttribute( adios2::IO &, std::string, + std::complex< long double > ) + { + throw std::runtime_error( + "[ADIOS2] Internal error: no support for long double complex attribute types" ); + } + + static void + readAttribute( adios2::IO &, std::string, + std::shared_ptr< Attribute::resource > ) + { + throw std::runtime_error( + "[ADIOS2] Internal error: no support for long double complex attribute types" ); + } + }; + + template< > struct AttributeTypes< std::vector< std::complex< long double > > > + { + using Attr = adios2::Attribute< std::complex< double > >; + using BasicType = double; + + static Attr createAttribute( adios2::IO &, std::string, + const std::vector< std::complex< long double > > & ) + { + throw std::runtime_error( + "[ADIOS2] Internal error: no support for long double complex vector attribute types" ); + } + + static void + readAttribute( adios2::IO &, std::string, + std::shared_ptr< Attribute::resource > ) + { + throw std::runtime_error( + "[ADIOS2] Internal error: no support for long double complex vector attribute types" ); + } + }; + template < typename T > struct AttributeTypes< std::vector< T > > { using Attr = adios2::Attribute< T >; @@ -504,6 +546,11 @@ namespace detail static constexpr bool validType = false; }; + // missing std::complex< long double > type in ADIOS2 v2.6.0 + template <> struct DatasetTypes< std::complex< long double > > + { + static constexpr bool validType = false; + }; template < typename T, size_t n > struct DatasetTypes< std::array< T, n > > { @@ -536,7 +583,7 @@ namespace detail /** * @brief Define a Variable of type T within the passed IO. - * + * * @param IO The adios2::IO object within which to define the * variable. The variable can later be retrieved from * the IO using the passed name. @@ -692,7 +739,7 @@ namespace detail private: /* - * ADIOS2 does not give direct access to its internal attribute and + * ADIOS2 does not give direct access to its internal attribute and * variable maps, but will instead give access to copies of them. * In order to avoid unnecessary copies, we buffer the returned map. * The downside of this is that we need to pay attention to invalidate @@ -702,7 +749,7 @@ namespace detail * been resolved * If false, the buffered map has been invalidated and needs to be * queried from ADIOS2 again. If true, the buffered map is equivalent to - * the map that would be returned by a call to + * the map that would be returned by a call to * IO::Available(Attributes|Variables). */ bool m_availableAttributesValid = false; diff --git a/include/openPMD/IO/HDF5/HDF5Auxiliary.hpp b/include/openPMD/IO/HDF5/HDF5Auxiliary.hpp index 26a5af8e25..abd1792afb 100644 --- a/include/openPMD/IO/HDF5/HDF5Auxiliary.hpp +++ b/include/openPMD/IO/HDF5/HDF5Auxiliary.hpp @@ -27,82 +27,103 @@ #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 +144,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 +241,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 94b44773da..84bbc4391e 100644 --- a/include/openPMD/IO/HDF5/HDF5IOHandlerImpl.hpp +++ b/include/openPMD/IO/HDF5/HDF5IOHandlerImpl.hpp @@ -66,7 +66,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/include/openPMD/IO/JSON/JSONIOHandlerImpl.hpp b/include/openPMD/IO/JSON/JSONIOHandlerImpl.hpp index 35ed7c8ef8..3367a35056 100644 --- a/include/openPMD/IO/JSON/JSONIOHandlerImpl.hpp +++ b/include/openPMD/IO/JSON/JSONIOHandlerImpl.hpp @@ -30,6 +30,7 @@ #include +#include #include #include #include @@ -145,6 +146,16 @@ namespace std return std::hash< shared_ptr< openPMD::File::FileState>> {}( s.fileState ); } }; + + // std::complex handling + template< class T > void to_json(nlohmann::json &j, const std::complex< T > &p) { + j = nlohmann::json {p.real(), p.imag()}; + } + + template< class T > void from_json(const nlohmann::json &j, std::complex< T > &p) { + p.real(j.at(0)); + p.imag(j.at(1)); + } } namespace openPMD diff --git a/include/openPMD/RecordComponent.hpp b/include/openPMD/RecordComponent.hpp index 80d8173866..299d9597af 100644 --- a/include/openPMD/RecordComponent.hpp +++ b/include/openPMD/RecordComponent.hpp @@ -273,7 +273,9 @@ RecordComponent::loadChunk(std::shared_ptr< T > data, Offset o, Extent e, double throw std::runtime_error("unitSI scaling during chunk loading not yet implemented"); Datatype dtype = determineDatatype(data); if( dtype != getDatatype() ) - if( !isSameInteger< T >( dtype ) && !isSameFloatingPoint< T >( dtype ) ) + if( !isSameInteger< T >( dtype ) && + !isSameFloatingPoint< T >( dtype ) && + !isSameComplexFloatingPoint< T >( dtype ) ) throw std::runtime_error("Type conversion during chunk loading not yet implemented"); uint8_t dim = getDimensionality(); diff --git a/include/openPMD/auxiliary/Memory.hpp b/include/openPMD/auxiliary/Memory.hpp index 1ce8c0a628..b54d30be74 100644 --- a/include/openPMD/auxiliary/Memory.hpp +++ b/include/openPMD/auxiliary/Memory.hpp @@ -23,6 +23,7 @@ #include "openPMD/Dataset.hpp" #include "openPMD/Datatype.hpp" +#include #include #include #include @@ -60,6 +61,21 @@ allocatePtr(Datatype dtype, uint64_t numPoints) data = new float[numPoints]; del = [](void* p){ delete[] static_cast< float* >(p); }; break; + case DT::VEC_CLONG_DOUBLE: + case DT::CLONG_DOUBLE: + data = new std::complex[numPoints]; + del = [](void* p){ delete[] static_cast< std::complex* >(p); }; + break; + case DT::VEC_CDOUBLE: + case DT::CDOUBLE: + data = new std::complex[numPoints]; + del = [](void* p){ delete[] static_cast< std::complex* >(p); }; + break; + case DT::VEC_CFLOAT: + case DT::CFLOAT: + data = new std::complex[numPoints]; + del = [](void* p){ delete[] static_cast< std::complex* >(p); }; + break; case DT::VEC_SHORT: case DT::SHORT: data = new short[numPoints]; diff --git a/include/openPMD/backend/Attribute.hpp b/include/openPMD/backend/Attribute.hpp index 25a2c1b9ec..a5effce539 100644 --- a/include/openPMD/backend/Attribute.hpp +++ b/include/openPMD/backend/Attribute.hpp @@ -25,6 +25,7 @@ #include #include +#include #include #include #include @@ -52,6 +53,7 @@ class Attribute : short, int, long, long long, unsigned short, unsigned int, unsigned long, unsigned long long, float, double, long double, + std::complex< float >, std::complex< double >, std::complex< long double >, std::string, std::vector< char >, std::vector< short >, @@ -66,6 +68,9 @@ class Attribute : std::vector< float >, std::vector< double >, std::vector< long double >, + std::vector< std::complex< float > >, + std::vector< std::complex< double > >, + std::vector< std::complex< long double > >, std::vector< std::string >, std::array< double, 7 >, bool > @@ -171,6 +176,12 @@ getCast( Attribute const & a ) return DoConvert{}(pvalue_d); else if(auto pvalue_ld = variantSrc::get_if< long double >( &v ) ) return DoConvert{}(pvalue_ld); + else if(auto pvalue_cf = variantSrc::get_if< std::complex< float > >( &v ) ) + return DoConvert, U>{}(pvalue_cf); + else if(auto pvalue_cd = variantSrc::get_if< std::complex< double > >( &v ) ) + return DoConvert, U>{}(pvalue_cd); + else if(auto pvalue_cld = variantSrc::get_if< std::complex< long double > >( &v ) ) + return DoConvert, U>{}(pvalue_cld); else if(auto pvalue_str = variantSrc::get_if< std::string >( &v ) ) return DoConvert{}(pvalue_str); // vector @@ -200,6 +211,12 @@ getCast( Attribute const & a ) return DoConvert, U>{}(pvalue_vd); else if(auto pvalue_vld = variantSrc::get_if< std::vector< long double > >( &v ) ) return DoConvert, U>{}(pvalue_vld); + else if(auto pvalue_vcf = variantSrc::get_if< std::vector< std::complex< float > > >( &v ) ) + return DoConvert >, U>{}(pvalue_vcf); + else if(auto pvalue_vcd = variantSrc::get_if< std::vector< std::complex< double > > >( &v ) ) + return DoConvert >, U>{}(pvalue_vcd); + else if(auto pvalue_vcld = variantSrc::get_if< std::vector< std::complex< long double > > >( &v ) ) + return DoConvert >, U>{}(pvalue_vcld); else if(auto pvalue_vstr = variantSrc::get_if< std::vector< std::string > >( &v ) ) return DoConvert, U>{}(pvalue_vstr); // extra diff --git a/include/openPMD/binding/python/Numpy.hpp b/include/openPMD/binding/python/Numpy.hpp index b3f29a3dc1..0f6bf9bc65 100644 --- a/include/openPMD/binding/python/Numpy.hpp +++ b/include/openPMD/binding/python/Numpy.hpp @@ -57,6 +57,12 @@ namespace openPMD return Datatype::ULONG; else if( dt.is(pybind11::dtype("ulonglong")) ) return Datatype::ULONGLONG; + else if( dt.is(pybind11::dtype("clongdouble")) ) + return Datatype::CLONG_DOUBLE; + else if( dt.is(pybind11::dtype("cdouble")) ) + return Datatype::CDOUBLE; + else if( dt.is(pybind11::dtype("csingle")) ) + return Datatype::CFLOAT; else if( dt.is(pybind11::dtype("longdouble")) ) return Datatype::LONG_DOUBLE; else if( dt.is(pybind11::dtype("double")) ) @@ -65,8 +71,10 @@ namespace openPMD return Datatype::FLOAT; else if( dt.is(pybind11::dtype("bool")) ) return Datatype::BOOL; - else + else { + pybind11::print(dt); throw std::runtime_error("Datatype '...' not known in 'dtype_from_numpy'!"); // _s.format(dt) + } } /** Return openPMD::Datatype from py::buffer_info::format @@ -79,7 +87,7 @@ namespace openPMD // refs: // https://docs.scipy.org/doc/numpy-1.15.0/reference/arrays.interface.html // https://docs.python.org/3/library/struct.html#format-characters - // std::cout << " scalar type '" << buf.format << "'" << std::endl; + // std::cout << " scalar type '" << fmt << "'" << std::endl; // typestring: encoding + type + number of bytes if( fmt.find("?") != std::string::npos ) return DT::BOOL; @@ -103,6 +111,12 @@ namespace openPMD return DT::ULONG; else if( fmt.find("Q") != std::string::npos ) return DT::ULONGLONG; + else if( fmt.find("Zf") != std::string::npos ) + return DT::CFLOAT; + else if( fmt.find("Zd") != std::string::npos ) + return DT::CDOUBLE; + else if( fmt.find("Zg") != std::string::npos ) + return DT::CLONG_DOUBLE; else if( fmt.find("f") != std::string::npos ) return DT::FLOAT; else if( fmt.find("d") != std::string::npos ) @@ -179,6 +193,18 @@ namespace openPMD case DT::VEC_LONG_DOUBLE: return pybind11::dtype("longdouble"); break; + case DT::CFLOAT: + case DT::VEC_CFLOAT: + return pybind11::dtype("csingle"); + break; + case DT::CDOUBLE: + case DT::VEC_CDOUBLE: + return pybind11::dtype("cdouble"); + break; + case DT::CLONG_DOUBLE: + case DT::VEC_CLONG_DOUBLE: + return pybind11::dtype("clongdouble"); + break; case DT::BOOL: return pybind11::dtype("bool"); // also "?" break; diff --git a/src/Datatype.cpp b/src/Datatype.cpp index 84faa7698d..ae3752263a 100644 --- a/src/Datatype.cpp +++ b/src/Datatype.cpp @@ -83,6 +83,15 @@ operator<<(std::ostream& os, openPMD::Datatype const & d) case DT::LONG_DOUBLE: os << "LONG_DOUBLE"; break; + case DT::CFLOAT: + os << "CFLOAT"; + break; + case DT::CDOUBLE: + os << "CDOUBLE"; + break; + case DT::CLONG_DOUBLE: + os << "CLONG_DOUBLE"; + break; case DT::STRING: os << "STRING"; break; @@ -125,6 +134,15 @@ operator<<(std::ostream& os, openPMD::Datatype const & d) case DT::VEC_LONG_DOUBLE: os << "VEC_LONG_DOUBLE"; break; + case DT::VEC_CFLOAT: + os << "VEC_CFLOAT"; + break; + case DT::VEC_CDOUBLE: + os << "VEC_CDOUBLE"; + break; + case DT::VEC_CLONG_DOUBLE: + os << "VEC_CLONG_DOUBLE"; + break; case DT::VEC_STRING: os << "VEC_STRING"; break; @@ -203,6 +221,18 @@ operator<<(std::ostream& os, openPMD::Datatype const & d) "LONG_DOUBLE", Datatype::LONG_DOUBLE }, + { + "CFLOAT", + Datatype::CFLOAT + }, + { + "CDOUBLE", + Datatype::CDOUBLE + }, + { + "CLONG_DOUBLE", + Datatype::CLONG_DOUBLE + }, { "STRING", Datatype::STRING @@ -259,6 +289,18 @@ operator<<(std::ostream& os, openPMD::Datatype const & d) "VEC_LONG_DOUBLE", Datatype::VEC_LONG_DOUBLE }, + { + "VEC_CFLOAT", + Datatype::VEC_CFLOAT + }, + { + "VEC_CDOUBLE", + Datatype::VEC_CDOUBLE + }, + { + "VEC_CLONG_DOUBLE", + Datatype::VEC_CLONG_DOUBLE + }, { "VEC_STRING", Datatype::VEC_STRING @@ -314,6 +356,9 @@ operator<<(std::ostream& os, openPMD::Datatype const & d) Datatype::FLOAT, Datatype::DOUBLE, Datatype::LONG_DOUBLE, + Datatype::CFLOAT, + Datatype::CDOUBLE, + Datatype::CLONG_DOUBLE, Datatype::STRING, Datatype::VEC_CHAR, Datatype::VEC_SHORT, @@ -328,6 +373,9 @@ operator<<(std::ostream& os, openPMD::Datatype const & d) Datatype::VEC_FLOAT, Datatype::VEC_DOUBLE, Datatype::VEC_LONG_DOUBLE, + Datatype::VEC_CFLOAT, + Datatype::VEC_CDOUBLE, + Datatype::VEC_CLONG_DOUBLE, Datatype::VEC_STRING, Datatype::ARR_DBL_7, Datatype::BOOL, @@ -363,7 +411,7 @@ operator<<(std::ostream& os, openPMD::Datatype const & d) if (it != map.end()) { return it->second; } else { - std::cerr << "Encountered non-basice type " << dt << ", aborting." + std::cerr << "Encountered non-basic type " << dt << ", aborting." << std::endl; throw std::runtime_error("toVectorType: passed non-basic type."); } diff --git a/src/IO/ADIOS/ADIOS2Auxiliary.cpp b/src/IO/ADIOS/ADIOS2Auxiliary.cpp index 8c1ff70577..e63f6cf59a 100644 --- a/src/IO/ADIOS/ADIOS2Auxiliary.cpp +++ b/src/IO/ADIOS/ADIOS2Auxiliary.cpp @@ -93,6 +93,9 @@ namespace detail { "float", Datatype::FLOAT }, { "double", Datatype::DOUBLE }, { "long double", Datatype::LONG_DOUBLE }, + { "float complex", Datatype::CFLOAT }, + { "double complex", Datatype::CDOUBLE }, + { "long double complex", Datatype::CLONG_DOUBLE }, // does not exist as of 2.6.0 but might come later { "uint8_t", Datatype::UCHAR }, { "int8_t", Datatype::CHAR }, { "uint16_t", determineDatatype< uint16_t >() }, diff --git a/src/IO/ADIOS/ADIOS2IOHandler.cpp b/src/IO/ADIOS/ADIOS2IOHandler.cpp index 909be4321e..3dee8a113d 100644 --- a/src/IO/ADIOS/ADIOS2IOHandler.cpp +++ b/src/IO/ADIOS/ADIOS2IOHandler.cpp @@ -771,8 +771,9 @@ namespace detail { } - template < typename T > - void DatasetReader::operator( )( detail::BufferedGet & bp, adios2::IO & IO, + template < typename T> + void + DatasetReader::operator( )( detail::BufferedGet & bp, adios2::IO & IO, adios2::Engine & engine, std::string const & fileName ) { diff --git a/src/IO/ADIOS/CommonADIOS1IOHandler.cpp b/src/IO/ADIOS/CommonADIOS1IOHandler.cpp index fdc7c13268..cb4f861279 100644 --- a/src/IO/ADIOS/CommonADIOS1IOHandler.cpp +++ b/src/IO/ADIOS/CommonADIOS1IOHandler.cpp @@ -19,6 +19,7 @@ * If not, see . */ #include +#include #include #include @@ -185,6 +186,23 @@ CommonADIOS1IOHandlerImpl::flush_attribute(int64_t group, std::string const& nam *ptr = att.get< long double >(); break; } + case DT::CFLOAT: + { + auto ptr = reinterpret_cast< std::complex< float >* >(values.get()); + *ptr = att.get< std::complex< float > >(); + break; + } + case DT::CDOUBLE: + { + auto ptr = reinterpret_cast< std::complex< double >* >(values.get()); + *ptr = att.get< std::complex< double > >(); + break; + } + case DT::CLONG_DOUBLE: + { + throw std::runtime_error("[ADIOS1] Unknown Attribute datatype (CLONG_DOUBLE)"); + break; + } case DT::STRING: { auto const & v = att.get< std::string >(); @@ -296,6 +314,16 @@ CommonADIOS1IOHandlerImpl::flush_attribute(int64_t group, std::string const& nam ptr[i] = vec[i]; break; } + /* not supported by ADIOS 1.13.1: + * https://github.com/ornladios/ADIOS/issues/212 + */ + case DT::VEC_CFLOAT: + case DT::VEC_CDOUBLE: + case DT::VEC_CLONG_DOUBLE: + { + throw std::runtime_error("[ADIOS1] Arrays of complex attributes are not supported"); + break; + } case DT::VEC_STRING: { auto ptr = reinterpret_cast< char** >(values.get()); @@ -735,11 +763,15 @@ CommonADIOS1IOHandlerImpl::openDataset(Writable* writable, case adios_long_double: dtype = DT::LONG_DOUBLE; break; + case adios_complex: + dtype = DT::CFLOAT; + break; + case adios_double_complex: + dtype = DT::CDOUBLE; + break; case adios_string: case adios_string_array: - case adios_complex: - case adios_double_complex: default: throw unsupported_data_error("[ADIOS1] Datatype not implemented for ADIOS dataset writing"); } @@ -900,6 +932,8 @@ CommonADIOS1IOHandlerImpl::readDataset(Writable* writable, using DT = Datatype; case DT::DOUBLE: case DT::FLOAT: + case DT::CDOUBLE: + case DT::CFLOAT: case DT::SHORT: case DT::INT: case DT::LONG: @@ -1013,14 +1047,18 @@ CommonADIOS1IOHandlerImpl::readAttribute(Writable* writable, case adios_long_double: size /= sizeof(long double); break; + case adios_complex: + size /= 8; + break; + case adios_double_complex: + size /= 16; + break; case adios_string: break; case adios_string_array: size /= sizeof(char*); break; - case adios_complex: - case adios_double_complex: default: throw unsupported_data_error( "[ADIOS1] readAttribute: Unsupported ADIOS1 attribute datatype '" + @@ -1198,6 +1236,14 @@ CommonADIOS1IOHandlerImpl::readAttribute(Writable* writable, dtype = DT::LONG_DOUBLE; a = Attribute(*reinterpret_cast< long double* >(data)); break; + case adios_complex: + dtype = DT::CFLOAT; + a = Attribute(*reinterpret_cast< std::complex* >(data)); + break; + case adios_double_complex: + dtype = DT::CDOUBLE; + a = Attribute(*reinterpret_cast< std::complex* >(data)); + break; case adios_string: { dtype = DT::STRING; @@ -1220,8 +1266,6 @@ CommonADIOS1IOHandlerImpl::readAttribute(Writable* writable, a = Attribute(vs); break; } - case adios_complex: - case adios_double_complex: default: throw unsupported_data_error( "[ADIOS1] readAttribute: Unsupported ADIOS1 attribute datatype '" + @@ -1372,6 +1416,9 @@ CommonADIOS1IOHandlerImpl::readAttribute(Writable* writable, a = Attribute(vld); break; } + /* not supported by ADIOS 1.13.1: VEC_CFLOAT, VEC_CDOUBLE, VEC_CLONG_DOUBLE + * https://github.com/ornladios/ADIOS/issues/212 + */ case adios_string: { dtype = DT::STRING; @@ -1394,8 +1441,6 @@ CommonADIOS1IOHandlerImpl::readAttribute(Writable* writable, break; } - case adios_complex: - case adios_double_complex: default: throw unsupported_data_error( "[ADIOS1] readAttribute: Unsupported ADIOS1 attribute datatype '" + diff --git a/src/IO/HDF5/HDF5IOHandler.cpp b/src/IO/HDF5/HDF5IOHandler.cpp index fef2ce4ecb..833c9960ea 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"; } } @@ -221,7 +250,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); @@ -245,7 +274,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; /* { @@ -257,11 +286,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; } @@ -269,9 +298,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, @@ -523,6 +558,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) ) @@ -756,6 +799,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; @@ -764,8 +814,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: @@ -822,11 +876,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 ) @@ -931,6 +987,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, @@ -1001,6 +1075,21 @@ HDF5IOHandlerImpl::writeAttribute(Writable* writable, dataType, att.get< std::vector< long double > >().data()); break; + case DT::VEC_CFLOAT: + status = H5Awrite(attribute_id, + dataType, + att.get< std::vector< std::complex< float > > >().data()); + break; + case DT::VEC_CDOUBLE: + status = H5Awrite(attribute_id, + dataType, + att.get< std::vector< std::complex< double > > >().data()); + break; + case DT::VEC_CLONG_DOUBLE: + status = H5Awrite(attribute_id, + dataType, + att.get< std::vector< std::complex< long double > > >().data()); + break; case DT::VEC_STRING: { auto vs = att.get< std::vector< std::string > >(); @@ -1032,11 +1121,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"); @@ -1085,8 +1171,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: @@ -1106,6 +1196,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, @@ -1303,6 +1399,17 @@ HDF5IOHandlerImpl::readAttribute(Writable* writable, throw unsupported_data_error("[HDF5] Unsupported attribute enumeration"); } else if( H5Tget_class(attr_type) == H5T_COMPOUND ) { + bool isComplexType = false; + if( H5Tget_nmembers(attr_type) == 2 ) + { + char* m0 = H5Tget_member_name(attr_type, 0); + char* m1 = H5Tget_member_name(attr_type, 1); + if( (strncmp("r" , m0, 4) == 0) && (strncmp("i", m1, 5) == 0) ) + isComplexType = true; + H5free_memory(m1); + H5free_memory(m0); + } + // re-implement legacy libSplash attributes for ColDim // see: include/splash/basetypes/ColTypeDim.hpp bool isLegacyLibSplashAttr = ( @@ -1327,6 +1434,29 @@ HDF5IOHandlerImpl::readAttribute(Writable* writable, attr_type, vc.data()); a = Attribute(vc); + } else if( isComplexType ) + { + size_t complexSize = H5Tget_member_offset(attr_type, 1); + if( complexSize == sizeof(float) ) + { + std::complex< float > cf; + status = H5Aread(attr_id, attr_type, &cf); + a = Attribute(cf); + } + else if( complexSize == sizeof(double) ) + { + std::complex< double > cd; + status = H5Aread(attr_id, attr_type, &cd); + a = Attribute(cd); + } + else if( complexSize == sizeof(long double) ) + { + std::complex< long double > cld; + status = H5Aread(attr_id, attr_type, &cld); + a = Attribute(cld); + } + else + throw unsupported_data_error("[HDF5] Unknow complex type representation"); } else throw unsupported_data_error("[HDF5] Compound attribute type not supported"); @@ -1439,6 +1569,27 @@ HDF5IOHandlerImpl::readAttribute(Writable* writable, attr_type, vld.data()); a = Attribute(vld); + } else if( H5Tequal(attr_type, m_H5T_CFLOAT) ) + { + std::vector< std::complex< float > > vcf(dims[0], 0); + status = H5Aread(attr_id, + attr_type, + vcf.data()); + a = Attribute(vcf); + } else if( H5Tequal(attr_type, m_H5T_CDOUBLE) ) + { + std::vector< std::complex< double > > vcd(dims[0], 0); + status = H5Aread(attr_id, + attr_type, + vcd.data()); + a = Attribute(vcd); + } else if( H5Tequal(attr_type, m_H5T_CLONG_DOUBLE) ) + { + std::vector< std::complex< long double > > vcld(dims[0], 0); + status = H5Aread(attr_id, + attr_type, + vcld.data()); + a = Attribute(vcld); } else if( H5Tget_class(attr_type) == H5T_STRING ) { std::vector< std::string > vs; diff --git a/src/IO/JSON/JSONIOHandlerImpl.cpp b/src/IO/JSON/JSONIOHandlerImpl.cpp index 2497ff9508..f45e7c6c08 100644 --- a/src/IO/JSON/JSONIOHandlerImpl.cpp +++ b/src/IO/JSON/JSONIOHandlerImpl.cpp @@ -1335,6 +1335,9 @@ namespace openPMD Datatype::FLOAT, Datatype::DOUBLE, Datatype::LONG_DOUBLE, + Datatype::CFLOAT, + Datatype::CDOUBLE, + Datatype::CLONG_DOUBLE, Datatype::BOOL }; for( auto it = std::begin( datatypes ); diff --git a/src/RecordComponent.cpp b/src/RecordComponent.cpp index 1ee776a555..17273c0131 100644 --- a/src/RecordComponent.cpp +++ b/src/RecordComponent.cpp @@ -27,6 +27,7 @@ #include #include #include +#include namespace openPMD @@ -188,6 +189,15 @@ RecordComponent::readBase() case DT::FLOAT: makeConstant(a.get< float >()); break; + case DT::CLONG_DOUBLE: + makeConstant(a.get< std::complex< long double > >()); + break; + case DT::CDOUBLE: + makeConstant(a.get< std::complex< double > >()); + break; + case DT::CFLOAT: + makeConstant(a.get< std::complex< float > >()); + break; case DT::SHORT: makeConstant(a.get< short >()); break; diff --git a/src/backend/Attributable.cpp b/src/backend/Attributable.cpp index 6403a42d1e..7eb402c1f2 100644 --- a/src/backend/Attributable.cpp +++ b/src/backend/Attributable.cpp @@ -23,6 +23,7 @@ #include #include +#include namespace openPMD @@ -219,6 +220,15 @@ Attributable::readAttributes() case DT::LONG_DOUBLE: setAttribute(att, a.get< long double >()); break; + case DT::CFLOAT: + setAttribute(att, a.get< std::complex< float > >()); + break; + case DT::CDOUBLE: + setAttribute(att, a.get< std::complex< double > >()); + break; + case DT::CLONG_DOUBLE: + setAttribute(att, a.get< std::complex< long double > >()); + break; case DT::STRING: setAttribute(att, a.get< std::string >()); break; @@ -261,6 +271,15 @@ Attributable::readAttributes() case DT::VEC_LONG_DOUBLE: setAttribute(att, a.get< std::vector< long double > >()); break; + case DT::VEC_CFLOAT: + setAttribute(att, a.get< std::vector< std::complex< float > > >()); + break; + case DT::VEC_CDOUBLE: + setAttribute(att, a.get< std::vector< std::complex< double > > >()); + break; + case DT::VEC_CLONG_DOUBLE: + setAttribute(att, a.get< std::vector< std::complex< long double > > >()); + break; case DT::VEC_STRING: setAttribute(att, a.get< std::vector< std::string > >()); break; diff --git a/src/binding/python/Attributable.cpp b/src/binding/python/Attributable.cpp index 1f1cdd9e5f..36af52bbce 100644 --- a/src/binding/python/Attributable.cpp +++ b/src/binding/python/Attributable.cpp @@ -27,9 +27,10 @@ #include "openPMD/auxiliary/Variant.hpp" #include "openPMD/binding/python/Numpy.hpp" +#include +#include #include #include -#include // std::variant @@ -113,6 +114,15 @@ bool setAttributeFromBufferInfo( case DT::LONG_DOUBLE: return attr.setAttribute( key, *static_cast(buf.ptr) ); break; + case DT::CFLOAT: + return attr.setAttribute( key, *static_cast*>(buf.ptr) ); + break; + case DT::CDOUBLE: + return attr.setAttribute( key, *static_cast*>(buf.ptr) ); + break; + case DT::CLONG_DOUBLE: + return attr.setAttribute( key, *static_cast*>(buf.ptr) ); + break; default: throw std::runtime_error("set_attribute: Unknown " "Python type '" + buf.format + @@ -155,6 +165,7 @@ bool setAttributeFromBufferInfo( static_cast(buf.ptr) + buf.size ) ); else */ + // std::cout << "+++++++++++ BUFFER: " << buf.format << std::endl; if( buf.format.find("b") != std::string::npos ) return attr.setAttribute( key, std::vector( @@ -215,6 +226,24 @@ bool setAttributeFromBufferInfo( static_cast(buf.ptr), static_cast(buf.ptr) + buf.size ) ); + else if( buf.format.find("Zf") != std::string::npos ) + return attr.setAttribute( key, + std::vector>( + static_cast*>(buf.ptr), + static_cast*>(buf.ptr) + buf.size + ) ); + else if( buf.format.find("Zd") != std::string::npos ) + return attr.setAttribute( key, + std::vector>( + static_cast*>(buf.ptr), + static_cast*>(buf.ptr) + buf.size + ) ); + else if( buf.format.find("Zg") != std::string::npos ) + return attr.setAttribute( key, + std::vector>( + static_cast*>(buf.ptr), + static_cast*>(buf.ptr) + buf.size + ) ); else if( buf.format.find("f") != std::string::npos ) return attr.setAttribute( key, std::vector( @@ -307,7 +336,7 @@ void init_Attributable(py::module &m) { // .def("set_attribute", &Attributable::setAttribute< std::vector< char > >) .def("set_attribute", &Attributable::setAttribute< std::vector< unsigned char > >) .def("set_attribute", &Attributable::setAttribute< std::vector< long > >) - .def("set_attribute", &Attributable::setAttribute< std::vector< double > >) + .def("set_attribute", &Attributable::setAttribute< std::vector< double > >) // TODO: this implicitly casts list of complex // probably affected by bug https://github.com/pybind/pybind11/issues/1258 .def("set_attribute", []( Attributable & attr, std::string const& key, std::vector< std::string > const& value ) { return attr.setAttribute( key, value ); diff --git a/src/binding/python/RecordComponent.cpp b/src/binding/python/RecordComponent.cpp index 4d9841ff98..d474f4ea3a 100644 --- a/src/binding/python/RecordComponent.cpp +++ b/src/binding/python/RecordComponent.cpp @@ -27,6 +27,7 @@ #include "openPMD/auxiliary/ShareRaw.hpp" #include "openPMD/binding/python/Numpy.hpp" +#include #include #include #include @@ -395,6 +396,12 @@ load_chunk(RecordComponent & r, Offset const & offset, Extent const & extent, st r.loadChunk(shareRaw((double*) a.mutable_data()), offset, extent); else if( r.getDatatype() == Datatype::FLOAT ) r.loadChunk(shareRaw((float*) a.mutable_data()), offset, extent); + else if( r.getDatatype() == Datatype::CLONG_DOUBLE ) + r.loadChunk>(shareRaw((std::complex*) a.mutable_data()), offset, extent); + else if( r.getDatatype() == Datatype::CDOUBLE ) + r.loadChunk>(shareRaw((std::complex*) a.mutable_data()), offset, extent); + else if( r.getDatatype() == Datatype::CFLOAT ) + r.loadChunk>(shareRaw((std::complex*) a.mutable_data()), offset, extent); else if( r.getDatatype() == Datatype::BOOL ) r.loadChunk(shareRaw((bool*) a.mutable_data()), offset, extent); else @@ -508,6 +515,15 @@ void init_RecordComponent(py::module &m) { case DT::LONG_DOUBLE: return rc.makeConstant( *static_cast(buf.ptr) ); break; + case DT::CFLOAT: + return rc.makeConstant( *static_cast*>(buf.ptr) ); + break; + case DT::CDOUBLE: + return rc.makeConstant( *static_cast*>(buf.ptr) ); + break; + case DT::CLONG_DOUBLE: + return rc.makeConstant( *static_cast*>(buf.ptr) ); + break; default: throw std::runtime_error("make_constant: " "Unknown Datatype!"); diff --git a/test/CoreTest.cpp b/test/CoreTest.cpp index 9e24efaf6a..5079cf6df4 100644 --- a/test/CoreTest.cpp +++ b/test/CoreTest.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -61,6 +62,12 @@ TEST_CASE( "attribute_dtype_test", "[core]" ) REQUIRE(Datatype::DOUBLE == a.dtype); a = Attribute(static_cast< long double >(0.)); REQUIRE(Datatype::LONG_DOUBLE == a.dtype); + a = Attribute(static_cast< std::complex< float > >(0.)); + REQUIRE(Datatype::CFLOAT == a.dtype); + a = Attribute(static_cast< std::complex< double > >(0.)); + REQUIRE(Datatype::CDOUBLE == a.dtype); + a = Attribute(static_cast< std::complex< long double > >(0.)); + REQUIRE(Datatype::CLONG_DOUBLE == a.dtype); a = Attribute(std::string("")); REQUIRE(Datatype::STRING == a.dtype); a = Attribute(std::vector< char >()); diff --git a/test/SerialIOTest.cpp b/test/SerialIOTest.cpp index 0a42b78352..d1e09557cb 100644 --- a/test/SerialIOTest.cpp +++ b/test/SerialIOTest.cpp @@ -13,6 +13,7 @@ #include #include #include +#include #include #include #include @@ -827,6 +828,82 @@ TEST_CASE( "write_test", "[serial]" ) } } +void test_complex(const std::string & backend) { + { + Series o = Series("../samples/serial_write_complex." + backend, AccessType::CREATE); + o.setAttribute("lifeIsComplex", std::complex(4.56, 7.89)); + o.setAttribute("butComplexFloats", std::complex(42.3, -99.3)); + if( backend != "bp" ) + o.setAttribute("longDoublesYouSay", std::complex(5.5, -4.55)); + + auto Cflt = o.iterations[0].meshes["Cflt"][RecordComponent::SCALAR]; + std::vector< std::complex > cfloats(3); + cfloats.at(0) = {1., 2.}; + cfloats.at(1) = {-3., 4.}; + cfloats.at(2) = {5., -6.}; + Cflt.resetDataset(Dataset(Datatype::CFLOAT, {cfloats.size()})); + Cflt.storeChunk(cfloats, {0}); + + auto Cdbl = o.iterations[0].meshes["Cdbl"][RecordComponent::SCALAR]; + std::vector< std::complex > cdoubles(3); + cdoubles.at(0) = {2., 1.}; + cdoubles.at(1) = {-4., 3.}; + cdoubles.at(2) = {6., -5.}; + Cdbl.resetDataset(Dataset(Datatype::CDOUBLE, {cdoubles.size()})); + Cdbl.storeChunk(cdoubles, {0}); + + std::vector< std::complex > cldoubles(3); + if( backend != "bp" ) + { + auto Cldbl = o.iterations[0].meshes["Cldbl"][RecordComponent::SCALAR]; + cldoubles.at(0) = {3., 2.}; + cldoubles.at(1) = {-5., 4.}; + cldoubles.at(2) = {7., -6.}; + Cldbl.resetDataset(Dataset(Datatype::CLONG_DOUBLE, {cldoubles.size()})); + Cldbl.storeChunk(cldoubles, {0}); + } + + o.flush(); + } + + //! @todo clarify that complex data is not N+1 data in JSON + if( backend != "json" ) + { + Series i = Series("../samples/serial_write_complex." + backend, AccessType::READ_ONLY); + REQUIRE(i.getAttribute("lifeIsComplex").get< std::complex >() == std::complex(4.56, 7.89)); + REQUIRE(i.getAttribute("butComplexFloats").get< std::complex >() == std::complex(42.3, -99.3)); + if( backend != "bp" ) { + REQUIRE(i.getAttribute("longDoublesYouSay").get >() == + std::complex(5.5, -4.55)); + } + + auto rcflt = i.iterations[0].meshes["Cflt"][RecordComponent::SCALAR].loadChunk< std::complex >(); + auto rcdbl = i.iterations[0].meshes["Cdbl"][RecordComponent::SCALAR].loadChunk< std::complex >(); + i.flush(); + + REQUIRE(rcflt.get()[1] == std::complex(-3., 4.)); + REQUIRE(rcdbl.get()[2] == std::complex(6, -5.)); + + if( backend != "bp" ) + { + auto rcldbl = i.iterations[0].meshes["Cldbl"][RecordComponent::SCALAR].loadChunk< std::complex >(); + i.flush(); + REQUIRE(rcldbl.get()[2] == std::complex(7., -6.)); + } + } +} + +TEST_CASE( "test_complex", "[serial]" ) +{ + // Notes: + // - ADIOS1 and ADIOS 2.6.0 have no complex long double + // - JSON read-back not distinguishable yet from N+1 shaped data set + for (auto const & t : getFileExtensions()) + { + test_complex(t); + } +} + inline void fileBased_add_EDpic(ParticleSpecies& e, uint64_t const num_particles) { diff --git a/test/python/unittest/API/APITest.py b/test/python/unittest/API/APITest.py index 2a0b2ad6a4..5b19a39025 100644 --- a/test/python/unittest/API/APITest.py +++ b/test/python/unittest/API/APITest.py @@ -158,6 +158,10 @@ def attributeRoundTrip(self, file_ending): series.set_attribute("single", np.single(1.234)) series.set_attribute("double", np.double(1.234567)) series.set_attribute("longdouble", np.longdouble(1.23456789)) + series.set_attribute("csingle", np.complex64(1.+2.j)) + series.set_attribute("cdouble", np.complex128(3.+4.j)) + if file_ending != "bp": + series.set_attribute("clongdouble", np.clongdouble(5.+6.j)) # array of ... series.set_attribute("arr_int16", (np.int16(23), np.int16(26), )) series.set_attribute("arr_int32", (np.int32(34), np.int32(37), )) @@ -183,6 +187,19 @@ def attributeRoundTrip(self, file_ending): series.set_attribute("l_double", [np.double(6.7), np.double(7.1)]) series.set_attribute("l_longdouble", [np.longdouble(7.8e9), np.longdouble(8.2e3)]) + # TODO: ComplexWarning: Casting complex values to real discards the + # imaginary part + # series.set_attribute("l_csingle", + # [np.csingle(5.6+7.8j), + # np.csingle(5.9+5.8j)]) + # series.set_attribute("l_cdouble", + # [np.complex_(6.7+6.8j), + # np.complex_(7.1+7.2j)]) + # if file_ending != "bp": + # series.set_attribute("l_clongdouble", + # [np.clongfloat(7.8e9-6.5e9j), + # np.clongfloat(8.2e3-9.1e3j)]) + # numpy.array of ... series.set_attribute("nparr_int16", np.array([234, 567], dtype=np.int16)) @@ -196,6 +213,21 @@ def attributeRoundTrip(self, file_ending): np.array([4.5, 6.7], dtype=np.double)) series.set_attribute("nparr_longdouble", np.array([8.9, 7.6], dtype=np.longdouble)) + # note: looks like ADIOS 1.13.1 cannot write arrays of complex + # as attributes (writes 1st value for single and crashes + # in write for complex double) + # https://github.com/ornladios/ADIOS/issues/212 + if series.backend != "ADIOS1": + series.set_attribute("nparr_csingle", + np.array([1.2 - 0.3j, 2.3 + 4.2j], + dtype=np.complex64)) + series.set_attribute("nparr_cdouble", + np.array([4.5 + 1.1j, 6.7 - 2.2j], + dtype=np.complex128)) + if file_ending != "bp": + series.set_attribute("nparr_clongdouble", + np.array([8.9 + 7.8j, 7.6 + 9.2j], + dtype=np.clongdouble)) # c_types # TODO remove the .value and handle types directly? @@ -239,6 +271,13 @@ def attributeRoundTrip(self, file_ending): 1.234567) self.assertAlmostEqual(series.get_attribute("longdouble"), 1.23456789) + np.testing.assert_almost_equal(series.get_attribute("csingle"), + np.complex64(1.+2.j)) + self.assertAlmostEqual(series.get_attribute("cdouble"), + 3.+4.j) + if file_ending != "bp": + self.assertAlmostEqual(series.get_attribute("clongdouble"), + 5.+6.j) # array of ... (will be returned as list) self.assertListEqual(series.get_attribute("arr_int16"), [np.int16(23), np.int16(26), ]) @@ -261,6 +300,14 @@ def attributeRoundTrip(self, file_ending): [np.double(6.7), np.double(7.1)]) self.assertListEqual(series.get_attribute("l_longdouble"), [np.longdouble(7.8e9), np.longdouble(8.2e3)]) + # TODO: l_csingle + # self.assertListEqual(series.get_attribute("l_cdouble"), + # [np.complex128(6.7 + 6.8j), + # np.double(7.1 + 7.2j)]) + # if file_ending != "bp": + # self.assertListEqual(series.get_attribute("l_clongdouble"), + # [np.clongdouble(7.8e9 - 6.5e9j), + # np.clongdouble(8.2e3 - 9.1e3j)]) # numpy.array of ... self.assertListEqual(series.get_attribute("nparr_int16"), @@ -275,6 +322,19 @@ def attributeRoundTrip(self, file_ending): series.get_attribute("nparr_double"), [4.5, 6.7]) np.testing.assert_almost_equal( series.get_attribute("nparr_longdouble"), [8.9, 7.6]) + # see https://github.com/ornladios/ADIOS/issues/212 + if series.backend != "ADIOS1": + np.testing.assert_almost_equal( + series.get_attribute("nparr_csingle"), + np.array([1.2 - 0.3j, 2.3 + 4.2j], + dtype=np.complex64)) + np.testing.assert_almost_equal( + series.get_attribute("nparr_cdouble"), + [4.5 + 1.1j, 6.7 - 2.2j]) + if file_ending != "bp": # not in ADIOS 1.13.1 nor ADIOS 2.6.0 + np.testing.assert_almost_equal( + series.get_attribute("nparr_clongdouble"), + [8.9 + 7.8j, 7.6 + 9.2j]) # TODO instead of returning lists, return all arrays as np.array? # self.assertEqual( # series.get_attribute("nparr_int16").dtype, np.int16) @@ -365,6 +425,20 @@ def makeConstantRoundTrip(self, file_ending): extent)) ms["longdouble"][SCALAR].make_constant(np.longdouble(1.23456789)) + ms["complex64"][SCALAR].reset_dataset( + DS(np.dtype("complex64"), extent)) + ms["complex64"][SCALAR].make_constant( + np.complex64(1.234 + 2.345j)) + ms["complex128"][SCALAR].reset_dataset( + DS(np.dtype("complex128"), extent)) + ms["complex128"][SCALAR].make_constant( + np.complex128(1.234567 + 2.345678j)) + if file_ending != "bp": + ms["clongdouble"][SCALAR].reset_dataset( + DS(np.dtype("clongdouble"), extent)) + ms["clongdouble"][SCALAR].make_constant( + np.clongdouble(1.23456789 + 2.34567890j)) + # flush and close file del series @@ -433,7 +507,16 @@ def makeConstantRoundTrip(self, file_ending): np.dtype('double')) self.assertTrue(ms["longdouble"][SCALAR].load_chunk(o, e).dtype == np.dtype('longdouble')) - + if file_ending != "json": + self.assertTrue(ms["complex64"][SCALAR].load_chunk(o, e).dtype + == np.dtype('complex64')) + self.assertTrue(ms["complex128"][SCALAR].load_chunk(o, e).dtype + == np.dtype('complex128')) + if file_ending != "bp": + self.assertTrue(ms["clongdouble"][SCALAR].load_chunk(o, e) + .dtype == np.dtype('clongdouble')) + + # FIXME: why does this even work w/o a flush() ? self.assertEqual(ms["int16"][SCALAR].load_chunk(o, e), np.int16(234)) self.assertEqual(ms["int32"][SCALAR].load_chunk(o, e), @@ -452,6 +535,13 @@ def makeConstantRoundTrip(self, file_ending): np.longdouble(1.23456789)) self.assertEqual(ms["double"][SCALAR].load_chunk(o, e), np.double(1.234567)) + self.assertEqual(ms["complex64"][SCALAR].load_chunk(o, e), + np.complex64(1.234 + 2.345j)) + self.assertEqual(ms["complex128"][SCALAR].load_chunk(o, e), + np.complex128(1.234567 + 2.345678j)) + if file_ending != "bp": + self.assertEqual(ms["clongdouble"][SCALAR].load_chunk(o, e), + np.clongdouble(1.23456789 + 2.34567890j)) def testConstantRecords(self): for ext in io.file_extensions: