How should pynwb support complex dtypes?

[pauladkisson]

Problem Description:
When processing time series data containing complex values (real and imaginary parts), Neurosift visualization failed (see flatironinstitute/neurosift#219). However, PyNWB successfully wrote the data, passed validation checks, and the NWB inspector showed no errors.

Technical Details:

Time series data contained complex data type (np.complex128)
PyNWB wrote, validated, and inspected the data without errors
Neurosift visualization engine failed when encountering the complex values

Decision Needed:
Should PyNWB:

Add official support of complex data types (np.complex128), or
Modify validation to throw an error when complex data types are detected?

[rly]

Complex dtype is technically not a supported dtype in NWB, even though those Python tools support writing and validating it. Here is a list of supported dtypes: https://schema-language.readthedocs.io/en/latest/description.html#dtype

Could you provide more information about the use case for storing processing time series data containing complex values? What do these data represent?

We could add official support for complex dtype, assuming that both the HDF5 and Zarr backends, and the Python, MATLAB APIs can handle it.

[pauladkisson]

Could you provide more information about the use case for storing processing time series data containing complex values? What do these data represent?

Unfortunately the case in which this came up turned out to be a bug, so I don't have a concrete, legitimate use case. I just thought it was worth flagging for discussion. @bendichter suggested that complex dtype data shouldn't pass validation so that these sorts of bugs are caught earlier.

[rly]

I agree that for now, we should raise an error on init and during validation for complex dtype data since it is not supported.

[bendichter]

Yes, and I'll add that if we want to support complex numbers in hdmf in the future the path is a bit more complicated. HDF5 does not natively support complex values (I don't know about Zarr but I doubt it). h5py has a way of mapping from the numpy complex data type to a custom data type that it can interpret as complex. This involves a single dataset where each element is complex with a length of 2, rather than having two datasets. If we were to support complex numbers in the schema language, I don't know if this would be the most convenient solution. It seems to me it would be much easier to just do something like:

groups:
    - complex:
        datasets:
            - real:
                dtype: float
            - imag:
                dtype: float

which would not require any changes to the schema language

[oruebel]

Essentially, I believe there are 3 main options to represent complex numbers using the schema language:

1. Encode the complex number in the dtype by using a compound type dtype=(float,float) such that the complex numbers are stored in a single dataset
2. Encode the complex number in the shape of the dataset by using an additional dimension of the dataset to represent the values. I.e., similar to how an RGB image is stored as shape=(x, y, 3), complex number could be stored in the last dimension shape=(...,2)
3. Store the components of the complex number as separate datasets as proposed above

It ultimately depends on the use case, which option is preferred:

• Option 1 (compound type) has the advantages that it represents the complex number as a single unit, keeping real and imag parts together, and making it easy to read, write, and compute on the values as a single complex number unit. Conversely, if the use case requires independent computations on the real and imag components then the compound type is not ideal, because the the two components have to be written and read as a unit
• Option 3 (two datasets) is in many ways the opposite of Option 1. This option is ideal when we need to write, read, and compute on the real and imag components independent of each other, but is inefficient and inconvenient when one needs to act on the complex number as a single unit, because it requires manual read and assembly of data from two separate datasets.
• Option 2 (encode in shape) is in some sense a compromise between Option 1 and 3. It allows reading of the real and compound types independent of one another by slicing just one component of the last dimension, but it is not as efficient in doing so as Option 3. It also allows reading of the real and imag part in one operation but still requires some manual assembly as the components are stored as separate indicies in the last dimensions, so it is not as convenient as Option 1 for acting on the complex numbers as a unit.

All of these three options should work right now to explicitly model complex numbers in the schema and extensions.

With that said, the question then is, whether to: a) allow dtype=complex64 and dtype=complex128 in the schema, b) how to store these dtypes, and c) whether to allow complex numbers to be used in the PyNWB and MatNWB APIs when the dtype is numeric or any.

• Re a): Allowing dtype=complex64 and dtype=complex128 as new dtypes in the schema would in principle fine, but I don't think this is critical at the moment to add support for this in the schema language given the amount of effort required and given that Options 1-3 above are all already available for users to explicitly create this.
• Re b): In my opinion, type=complex128 in terms of schema would most naturally map to Option 1 above since it: 1) states to encode the values in the dtype and 2) having the dtype modify the shape or even create a new dataset would lead to complex side-effects that would be very complex to implement and that would likely lead to collisions in the schema. I.e., I think dtype=complex128 would essentially just be a shorthand for

dtype:
  - doc: The real part of the complex128 number
    dtype: float64
    name: real
  - doc: The imaginary part of the complex128 number
    dtype: float64
    name: imag

• Re c): Because complex numbers are not natively supported by either HDF5 nor Zarr, I think the effort required to implicitly map compound data types is likely quite high. This works right now, mainly because h5py implements its own custom way of mapping complex numbers to a non-standard dtype. But that process is non-standard. Options 3 (2 separate datasets) is impractical to implement when doing an implicit mapping because it could create collisions with other datasets that are already specified in the schema. Option 2 (encoding in the shape), would require significant custom logic in the I/O backend for both read and write. Option 1 (use a compound dtype) is essentially what h5py already does. At the end of the day, all of these options for implicitly mapping complex numbers leave us with custom logic that all APIs would need to implement and, hence, have a high burden for code maintenance and high risk for creating issues between APIs.

In my opinion, we already have options for users to explicitly model complex valued data, since Options 1,2, and 3 are already possible to do in a schema. With this, it is possible to explicitly model complex numbers in extensions. At the same time, the effort to support np.complex128 to implicitly store complex number in TimeSeries seems very high. Because of this I would recommend to:

• Not allow np.complex numbers as dtype for TimeSeries on write
• Create a short tutorial discussing how complex numbers can be implemented in extensions via the options above. In this way, users can create extensions if they need complex numbers.
• If complex numbers are common needed, then we could add a base class ComplexTimeSeries in the schema that would use one of the 3 options described above to explicitly model complex numbers, rather than creating API magic to implicitly map complex numbers

[bendichter]

Well put Oliver! I suppose a disadvantage of approach 3 is it does not enforce equal sized datasets so I am leaning towards 1 or 2 as a best practice.

Let's wait until we have an actual use case. The current situation was the result of a bug.

Community: comment here if you have a need for complex data!

[pauladkisson]

For a concrete minimal example (with mock data), see my example code below

from pynwb.testing.mock.base import mock_TimeSeries
from pynwb.testing.mock.file import mock_NWBFile
from pynwb import NWBHDF5IO
import numpy as np

def main():
    nwbfile = mock_NWBFile()
    complex_data = np.array([1 + 1j, 2 + 2j, 3 + 3j], dtype=np.complex128)
    time_series = mock_TimeSeries(data=complex_data, nwbfile=nwbfile)

    with NWBHDF5IO('complex_data.nwb', 'w') as io:
        io.write(nwbfile)
    
    with NWBHDF5IO('complex_data.nwb', 'r') as io:
        nwbfile = io.read()
        print(nwbfile.acquisition["TimeSeries"].data[:]) # [1.+1.j 2.+2.j 3.+3.j]
        print(nwbfile.acquisition["TimeSeries"].data.dtype) # complex128

if __name__ == "__main__":
    main()

When running pynwb-validate I get the following:

pynwb-validate complex_data.nwb
Validating complex_data.nwb against cached namespace information using namespace 'core'.
 - no errors found.

And from nwbinspector,

nwbinspector complex_data.nwb                                                                                              
/usr/local/lib/python3.9/site-packages/hdmf/utils.py:668: UserWarning: Ignoring cached namespace 'core' version 2.8.0 because version 2.5.0 is already loaded.
  return func(args[0], **pargs)
                                                                                                                                              

**************************************************
NWBInspector Report Summary

Timestamp: 2025-05-19 14:25:01.273179-07:00
Platform: macOS-13.4-x86_64-i386-64bit
NWBInspector version: 0.4.29

Found 7 issues over 1 files:
       7 - BEST_PRACTICE_SUGGESTION
**************************************************


0  BEST_PRACTICE_SUGGESTION
===========================

0.0  complex_data.nwb: check_description - 'TimeSeries' object at location '/acquisition/TimeSeries'
       Message: Description ('no description') is a placeholder.

0.1  complex_data.nwb: check_session_start_time_old_date - 'NWBFile' object at location '/'
       Message: The session_start_time (1970-01-01 00:00:00+00:00) may not be set to the true date of the recording.

0.2  complex_data.nwb: check_experimenter_exists - 'NWBFile' object at location '/'
       Message: Experimenter is missing.

0.3  complex_data.nwb: check_experiment_description - 'NWBFile' object at location '/'
       Message: Experiment description is missing.

0.4  complex_data.nwb: check_institution - 'NWBFile' object at location '/'
       Message: Metadata /general/institution is missing.

0.5  complex_data.nwb: check_keywords - 'NWBFile' object at location '/'
       Message: Metadata /general/keywords is missing.

0.6  complex_data.nwb: check_subject_exists - 'NWBFile' object at location '/'
       Message: Subject is missing.

[bendichter]

@rly should we move this to the hdmf-dev org?