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: ornladios/ADIOS#212)
- [x] ADIOS2: (complex attributes: ornladios/ADIOS2#1908);  (no long double complex ornladios/ADIOS2#1907)
- [x] Python bindings

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:
 

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 <https://github.com/ornladios/ADIOS/issues/212>`_.
+
 
 Selected References
 -------------------

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 <map>
 #include <stdexcept>
 #include <string>
+#include <complex>
+
 
 
 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 > >(

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:

include/openPMD/IO/ADIOS/ADIOS2Auxiliary.hpp

@@ -25,6 +25,8 @@
 #if openPMD_HAVE_ADIOS2
 #include "openPMD/Datatype.hpp"
 #include <adios2.h>
+#include <complex>
+#include <stdexcept>
 #include <utility>
 #include <vector>
 
@@ -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 > >
     {