test_autobad.py

import subprocess
import os
import mne

mne.set_log_level(False)

fname = 'bids/sub-{subject:02d}/ses-{session:02d}/meg/sub-{subject:02d}_ses-{session:02d}_task-naming_run-{run:02d}_meg.fif'

for subject in range(1, 21):
    for session in range(1, 3):
        for run in range(1, 4):
            raw_fname = fname.format(subject=subject, session=session, run=run)
            print(f'Processing: {raw_fname}', flush=True)

            try:
                os.remove('/tmp/out.fif')
            except FileNotFoundError:
                pass

            output = subprocess.check_output([
                'maxfilter', '-v',
                '-autobad', 'on',
                '-badlimit', '7',
                '-ctc', 'off',
                '-cal', 'off',
                '-regularize', 'off',
                '-f', raw_fname,
                '-o', '/tmp/out.fif',
            ], stderr=subprocess.STDOUT)
            os.remove('/tmp/out.fif')

            # Parse the output to look for these lines:
            #  Static bad channels (2): 511 205
            #  Detected 2 flat channels: 2233 2031
            maxwell_bads = set()
            maxwell_flats = set()
            for line in output.split(b'\n'):
                line = line.strip()
                if line.startswith(b'Static bad channels'):
                    bads = line.decode('utf8').split(': ')[1].split(' ')
                    for bad in bads:
                        maxwell_bads.add(bad)
                elif line.startswith(b'Detected'):
                    flats = line.decode('utf8').split(': ')[1].split(' ')
                    for flat in flats:
                        maxwell_flats.add(flat)

            # MNE version
            raw = mne.io.read_raw_fif(raw_fname, preload=True)
            mne.channels.fix_mag_coil_types(raw.info)
            raw.resample(98).filter(0.1, None)

            mne.preprocessing.mark_flat(raw, bad_percent=0.01)
            mne_flats = set(raw.info['bads'])
            # Use flat channels detected by maxfilter. This script is about
            # testing maxwell_autobad, not mark_flat.
            #raw.info['bads'] = ['MEG%04d' % int(f) for f in maxwell_flats]

            # Make the algorithm detect bads/flats by itself
            raw.info['bads'] = []

            mne_bads = set(mne.preprocessing.maxwell_autobad(
                raw,
                origin=(0., 0., 0.04),
                regularize=None,
                bad_condition='ignore',
            ))

            print('maxwell flats:', maxwell_flats)
            print('mne flats:', mne_flats)
            print('maxwell bads:', maxwell_bads)
            print('mne bads:', mne_bads)
            print(flush=True)