MRG, ENH: Export evokeds MFF (#9406)

* ENH: Add export evokeds to MFF function

* ENH: Add export evokeds function

* Add tests for export evokeds

* DOC: Add export_evokeds to python reference

* DOC: Fix history argument type

* STY: Use logger for user message

* FIX: Nest pytz import

* Clarify MFF specific arguments

* DOC: Use shared warning message

* FIX: Test digitization data from categories.xml

* STY: Use regex match for pytest.raises

* STY: Use zipped lists

* Refactor assert statements

* ENH: Store MFF device type in info

* Determine device type from info

* Use kwargs for format-specific exports

* Use current time for record time

Because averaged files are absolute time agnostic,
we can simplify by using plugging in the current time
for the record time. This is also the behavior when
running an averaging tool in EGI Net Station software.

* Add MFF export funtion to python reference

* ENH: Reorganize

* FIX: Missed a file

* ENH: For newer mffpy

Co-authored-by: Eric Larson <[email protected]>

doc/export.rst

@@ -14,4 +14,6 @@ Exporting
    :toctree: generated/
 
    export_epochs
+   export_evokeds
+   export_evokeds_mff
    export_raw

mne/__init__.py

@@ -77,7 +77,8 @@
                           events_from_annotations)
 from .epochs import (BaseEpochs, Epochs, EpochsArray, read_epochs,
                      concatenate_epochs, make_fixed_length_epochs)
-from .evoked import Evoked, EvokedArray, read_evokeds, write_evokeds, combine_evoked
+from .evoked import (Evoked, EvokedArray, read_evokeds, write_evokeds,
+                     combine_evoked)
 from .label import (read_label, label_sign_flip,
                     write_label, stc_to_label, grow_labels, Label, split_label,
                     BiHemiLabel, read_labels_from_annot, write_labels_to_annot,

mne/epochs.py

@@ -1823,7 +1823,7 @@ def export(self, fname, fmt='auto', verbose=None):
         """Export Epochs to external formats.
 
         Supported formats: EEGLAB (set, uses :mod:`eeglabio`)
-        %(export_warning)s
+        %(export_warning)s :meth:`save` instead.
 
         Parameters
         ----------

mne/export/__init__.py

@@ -1 +1,2 @@
-from ._export import export_raw, export_epochs
+from ._export import export_raw, export_epochs, export_evokeds
+from ._egimff import export_evokeds_mff

mne/export/_egimff.py

@@ -0,0 +1,150 @@
+# -*- coding: utf-8 -*-
+# Authors: MNE Developers
+#
+# License: BSD (3-clause)
+
+import datetime
+
+import numpy as np
+
+from ..io.egi.egimff import _import_mffpy
+from ..io.pick import pick_types, pick_channels
+from ..utils import verbose
+
+
+@verbose
+def export_evokeds_mff(fname, evoked, history=None, *, verbose=None):
+    """Export evoked dataset to MFF.
+
+    Parameters
+    ----------
+    %(export_params_fname)s
+    evoked : list of Evoked instances
+        List of evoked datasets to export to one file. Note that the
+        measurement info from the first evoked instance is used, so be sure
+        that information matches.
+    history : None (default) | list of dict
+        Optional list of history entries (dictionaries) to be written to
+        history.xml. This must adhere to the format described in
+        mffpy.xml_files.History.content. If None, no history.xml will be
+        written.
+    %(verbose)s
+
+    Notes
+    -----
+    .. versionadded:: 0.24
+
+    Only EEG channels are written to the output file.
+    ``info['device_info']['type']`` must be a valid MFF recording device
+    (e.g. 'HydroCel GSN 256 1.0'). This field is automatically populated when
+    using MFF read functions.
+    """
+    mffpy = _import_mffpy('Export evokeds to MFF.')
+    import pytz
+    info = evoked[0].info
+    if np.round(info['sfreq']) != info['sfreq']:
+        raise ValueError('Sampling frequency must be a whole number. '
+                         f'sfreq: {info["sfreq"]}')
+    sampling_rate = int(info['sfreq'])
+
+    # Initialize writer
+    writer = mffpy.Writer(fname)
+    current_time = pytz.utc.localize(datetime.datetime.utcnow())
+    writer.addxml('fileInfo', recordTime=current_time)
+    try:
+        device = info['device_info']['type']
+    except (TypeError, KeyError):
+        raise ValueError('No device type. Cannot determine sensor layout.')
+    writer.add_coordinates_and_sensor_layout(device)
+
+    # Add EEG data
+    eeg_channels = pick_types(info, eeg=True, exclude=[])
+    eeg_bin = mffpy.bin_writer.BinWriter(sampling_rate)
+    for ave in evoked:
+        # Signals are converted to µV
+        block = (ave.data[eeg_channels] * 1e6).astype(np.float32)
+        eeg_bin.add_block(block, offset_us=0)
+    writer.addbin(eeg_bin)
+
+    # Add categories
+    categories_content = _categories_content_from_evokeds(evoked)
+    writer.addxml('categories', categories=categories_content)
+
+    # Add history
+    if history:
+        writer.addxml('historyEntries', entries=history)
+
+    writer.write()
+
+
+def _categories_content_from_evokeds(evoked):
+    """Return categories.xml content for evoked dataset."""
+    content = dict()
+    begin_time = 0
+    for ave in evoked:
+        # Times are converted to microseconds
+        sfreq = ave.info['sfreq']
+        duration = np.round(len(ave.times) / sfreq * 1e6).astype(int)
+        end_time = begin_time + duration
+        event_time = begin_time - np.round(ave.tmin * 1e6).astype(int)
+        eeg_bads = _get_bad_eeg_channels(ave.info)
+        content[ave.comment] = [
+            _build_segment_content(begin_time, end_time, event_time, eeg_bads,
+                                   name='Average', nsegs=ave.nave)
+        ]
+        begin_time += duration
+    return content
+
+
+def _get_bad_eeg_channels(info):
+    """Return a list of bad EEG channels formatted for categories.xml.
+
+    Given a list of only the EEG channels in file, return the indices of this
+    list (starting at 1) that correspond to bad channels.
+    """
+    if len(info['bads']) == 0:
+        return []
+    eeg_channels = pick_types(info, eeg=True, exclude=[])
+    bad_channels = pick_channels(info['ch_names'], info['bads'])
+    bads_elementwise = np.isin(eeg_channels, bad_channels)
+    return list(np.flatnonzero(bads_elementwise) + 1)
+
+
+def _build_segment_content(begin_time, end_time, event_time, eeg_bads,
+                           status='unedited', name=None, pns_bads=None,
+                           nsegs=None):
+    """Build content for a single segment in categories.xml.
+
+    Segments are sorted into categories in categories.xml. In a segmented MFF
+    each category can contain multiple segments, but in an averaged MFF each
+    category only contains one segment (the average).
+    """
+    channel_status = [{
+        'signalBin': 1,
+        'exclusion': 'badChannels',
+        'channels': eeg_bads
+    }]
+    if pns_bads:
+        channel_status.append({
+            'signalBin': 2,
+            'exclusion': 'badChannels',
+            'channels': pns_bads
+        })
+    content = {
+        'status': status,
+        'beginTime': begin_time,
+        'endTime': end_time,
+        'evtBegin': event_time,
+        'evtEnd': event_time,
+        'channelStatus': channel_status,
+    }
+    if name:
+        content['name'] = name
+    if nsegs:
+        content['keys'] = {
+            '#seg': {
+                'type': 'long',
+                'data': nsegs
+            }
+        }
+    return content

mne/export/_export.py

@@ -5,7 +5,8 @@
 
 import os.path as op
 
-from ..utils import verbose, _validate_type
+from ._egimff import export_evokeds_mff
+from ..utils import verbose, logger, _validate_type
 
 
 @verbose
@@ -78,6 +79,63 @@ def export_epochs(fname, epochs, fmt='auto', verbose=None):
         raise NotImplementedError('Export to BrainVision not implemented.')
 
 
+@verbose
+def export_evokeds(fname, evoked, fmt='auto', verbose=None):
+    """Export evoked dataset to external formats.
+
+    This function is a wrapper for format-specific export functions. The export
+    function is selected based on the inferred file format. For additional
+    options, use the format-specific functions.
+
+    Supported formats
+        MFF (mff, uses :func:`mne.export.export_evokeds_mff`)
+    %(export_warning)s :func:`mne.write_evokeds` instead.
+
+    Parameters
+    ----------
+    %(export_params_fname)s
+    evoked : Evoked instance, or list of Evoked instances
+        The evoked dataset, or list of evoked datasets, to export to one file.
+        Note that the measurement info from the first evoked instance is used,
+        so be sure that information matches.
+    fmt : 'auto' | 'mff'
+        Format of the export. Defaults to ``'auto'``, which will infer the
+        format from the filename extension. See supported formats above for
+        more information.
+    %(verbose)s
+
+    See Also
+    --------
+    mne.write_evokeds
+    mne.export.export_evokeds_mff
+
+    Notes
+    -----
+    .. versionadded:: 0.24
+    """
+    supported_export_formats = {
+        'mff': ('mff',),
+        'eeglab': ('set',),
+        'edf': ('edf',),
+        'brainvision': ('eeg', 'vmrk', 'vhdr',)
+    }
+    fmt = _infer_check_export_fmt(fmt, fname, supported_export_formats)
+
+    if not isinstance(evoked, list):
+        evoked = [evoked]
+
+    logger.info(f'Exporting evoked dataset to {fname}...')
+
+    if fmt == 'mff':
+        export_evokeds_mff(fname, evoked)
+    elif fmt == 'eeglab':
+        raise NotImplementedError('Export to EEGLAB not implemented.')
+    elif fmt == 'edf':
+        raise NotImplementedError('Export to EDF not implemented.')
+    elif fmt == 'brainvision':
+        raise NotImplementedError('Export to BrainVision not implemented.')
+
+
 def _infer_check_export_fmt(fmt, fname, supported_formats):
     """Infer export format from filename extension if auto.
 

mne/export/tests/test_export.py

@@ -11,10 +11,18 @@
 import numpy as np
 from numpy.testing import assert_allclose, assert_array_equal
 
-from mne import read_epochs_eeglab, Epochs
-from mne.tests.test_epochs import _get_data
+from mne import read_epochs_eeglab, Epochs, read_evokeds, read_evokeds_mff
+from mne.datasets import testing
+from mne.export import export_evokeds, export_evokeds_mff
 from mne.io import read_raw_fif, read_raw_eeglab
-from mne.utils import _check_eeglabio_installed
+from mne.utils import _check_eeglabio_installed, requires_version, object_diff
+from mne.tests.test_epochs import _get_data
+
+base_dir = op.join(op.dirname(__file__), '..', '..', 'io', 'tests', 'data')
+fname_evoked = op.join(base_dir, 'test-ave.fif')
+
+data_path = testing.data_path(download=False)
+egi_evoked_fname = op.join(data_path, 'EGI', 'test_egi_evoked.mff')
 
 
 @pytest.mark.skipif(not _check_eeglabio_installed(strict=False),
@@ -62,3 +70,75 @@ def test_export_epochs_eeglab(tmpdir, preload):
     assert epochs.event_id.keys() == epochs_read.event_id.keys()  # just keys
     assert_allclose(epochs.times, epochs_read.times)
     assert_allclose(epochs.get_data(), epochs_read.get_data())
+
+
+@requires_version('mffpy', '0.5.7')
+@testing.requires_testing_data
+@pytest.mark.parametrize('fmt', ('auto', 'mff'))
+@pytest.mark.parametrize('do_history', (True, False))
+def test_export_evokeds_to_mff(tmpdir, fmt, do_history):
+    """Test exporting evoked dataset to MFF."""
+    evoked = read_evokeds_mff(egi_evoked_fname)
+    export_fname = op.join(str(tmpdir), 'evoked.mff')
+    history = [
+        {
+            'name': 'Test Segmentation',
+            'method': 'Segmentation',
+            'settings': ['Setting 1', 'Setting 2'],
+            'results': ['Result 1', 'Result 2']
+        },
+        {
+            'name': 'Test Averaging',
+            'method': 'Averaging',
+            'settings': ['Setting 1', 'Setting 2'],
+            'results': ['Result 1', 'Result 2']
+        }
+    ]
+    if do_history:
+        export_evokeds_mff(export_fname, evoked, history=history)
+    else:
+        export_evokeds(export_fname, evoked)
+    # Drop non-EEG channels
+    evoked = [ave.drop_channels(['ECG', 'EMG']) for ave in evoked]
+    evoked_exported = read_evokeds_mff(export_fname)
+    assert len(evoked) == len(evoked_exported)
+    for ave, ave_exported in zip(evoked, evoked_exported):
+        # Compare infos
+        assert object_diff(ave_exported.info, ave.info) == ''
+        # Compare data
+        assert_allclose(ave_exported.data, ave.data)
+        # Compare properties
+        assert ave_exported.nave == ave.nave
+        assert ave_exported.kind == ave.kind
+        assert ave_exported.comment == ave.comment
+        assert_allclose(ave_exported.times, ave.times)
+
+
+@requires_version('mffpy', '0.5.7')
+@testing.requires_testing_data
+def test_export_to_mff_no_device():
+    """Test no device type throws ValueError."""
+    evoked = read_evokeds_mff(egi_evoked_fname, condition='Category 1')
+    evoked.info['device_info'] = None
+    with pytest.raises(ValueError, match='No device type.'):
+        export_evokeds('output.mff', evoked)
+
+
+@requires_version('mffpy', '0.5.7')
+def test_export_to_mff_incompatible_sfreq():
+    """Test non-whole number sampling frequency throws ValueError."""
+    evoked = read_evokeds(fname_evoked)
+    with pytest.raises(ValueError, match=f'sfreq: {evoked[0].info["sfreq"]}'):
+        export_evokeds('output.mff', evoked)
+
+
+@pytest.mark.parametrize('fmt,ext', [
+    ('EEGLAB', 'set'),
+    ('EDF', 'edf'),
+    ('BrainVision', 'eeg')
+])
+def test_export_evokeds_unsupported_format(fmt, ext):
+    """Test exporting evoked dataset to non-supported formats."""
+    evoked = read_evokeds(fname_evoked)
+    with pytest.raises(NotImplementedError, match=f'Export to {fmt} not imp'):
+        export_evokeds(f'output.{ext}', evoked)

mne/io/base.py

@@ -1455,7 +1455,7 @@ def export(self, fname, fmt='auto', verbose=None):
         """Export Raw to external formats.
 
         Supported formats: EEGLAB (set, uses :mod:`eeglabio`)
-        %(export_warning)s
+        %(export_warning)s :meth:`save` instead.
 
         Parameters
         ----------