Python proc_history._get_sss_rank函数代码示例

def test_maxfilter_get_rank():
    """test maxfilter rank lookup"""
    raw = io.read_raw_fif(raw_fname)
    mf = raw.info['proc_history'][0]['max_info']
    rank1 = mf['sss_info']['nfree']
    rank2 = _get_sss_rank(mf)
    assert_equal(rank1, rank2)

def test_rank_estimation():
    """Test raw rank estimation
    """
    iter_tests = itt.product([fif_fname, hp_fif_fname], ["norm", dict(mag=1e11, grad=1e9, eeg=1e5)])  # sss
    for fname, scalings in iter_tests:
        raw = Raw(fname)
        (_, picks_meg), (_, picks_eeg) = _picks_by_type(raw.info, meg_combined=True)
        n_meg = len(picks_meg)
        n_eeg = len(picks_eeg)

        raw = Raw(fname, preload=True)
        if "proc_history" not in raw.info:
            expected_rank = n_meg + n_eeg
        else:
            mf = raw.info["proc_history"][0]["max_info"]
            expected_rank = _get_sss_rank(mf) + n_eeg
        assert_array_equal(raw.estimate_rank(scalings=scalings), expected_rank)

        assert_array_equal(raw.estimate_rank(picks=picks_eeg, scalings=scalings), n_eeg)

        raw = Raw(fname, preload=False)
        if "sss" in fname:
            tstart, tstop = 0.0, 30.0
            raw.add_proj(compute_proj_raw(raw))
            raw.apply_proj()
        else:
            tstart, tstop = 10.0, 20.0

        raw.apply_proj()
        n_proj = len(raw.info["projs"])

        assert_array_equal(
            raw.estimate_rank(tstart=tstart, tstop=tstop, scalings=scalings),
            expected_rank - (1 if "sss" in fname else n_proj),
        )

def test_basic():
    """Test Maxwell filter basic version"""
    # Load testing data (raw, SSS std origin, SSS non-standard origin)
    with warnings.catch_warnings(record=True):  # maxshield
        raw = Raw(raw_fname, allow_maxshield=True).crop(0., 1., False)
        raw_err = Raw(raw_fname, proj=True, allow_maxshield=True)
        raw_erm = Raw(erm_fname, allow_maxshield=True)
    assert_raises(RuntimeError, maxwell_filter, raw_err)
    assert_raises(TypeError, maxwell_filter, 1.)  # not a raw
    assert_raises(ValueError, maxwell_filter, raw, int_order=20)  # too many

    n_int_bases = int_order ** 2 + 2 * int_order
    n_ext_bases = ext_order ** 2 + 2 * ext_order
    nbases = n_int_bases + n_ext_bases

    # Check number of bases computed correctly
    assert_equal(_get_n_moments([int_order, ext_order]).sum(), nbases)

    # Test SSS computation at the standard head origin
    raw_sss = maxwell_filter(raw, origin=mf_head_origin, regularize=None,
                             bad_condition='ignore')
    assert_meg_snr(raw_sss, Raw(sss_std_fname), 200., 1000.)
    py_cal = raw_sss.info['proc_history'][0]['max_info']['sss_cal']
    assert_equal(len(py_cal), 0)
    py_ctc = raw_sss.info['proc_history'][0]['max_info']['sss_ctc']
    assert_equal(len(py_ctc), 0)
    py_st = raw_sss.info['proc_history'][0]['max_info']['max_st']
    assert_equal(len(py_st), 0)
    assert_raises(RuntimeError, maxwell_filter, raw_sss)

    # Test SSS computation at non-standard head origin
    raw_sss = maxwell_filter(raw, origin=[0., 0.02, 0.02], regularize=None,
                             bad_condition='ignore')
    assert_meg_snr(raw_sss, Raw(sss_nonstd_fname), 250., 700.)

    # Test SSS computation at device origin
    sss_erm_std = Raw(sss_erm_std_fname)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg',
                             origin=mf_meg_origin, regularize=None,
                             bad_condition='ignore')
    assert_meg_snr(raw_sss, sss_erm_std, 100., 900.)
    for key in ('job', 'frame'):
        vals = [x.info['proc_history'][0]['max_info']['sss_info'][key]
                for x in [raw_sss, sss_erm_std]]
        assert_equal(vals[0], vals[1])

    # Check against SSS functions from proc_history
    sss_info = raw_sss.info['proc_history'][0]['max_info']
    assert_equal(_get_n_moments(int_order),
                 proc_history._get_sss_rank(sss_info))

    # Degenerate cases
    raw_bad = raw.copy()
    raw_bad.comp = True
    assert_raises(RuntimeError, maxwell_filter, raw_bad)
    del raw_bad
    assert_raises(ValueError, maxwell_filter, raw, coord_frame='foo')
    assert_raises(ValueError, maxwell_filter, raw, origin='foo')
    assert_raises(ValueError, maxwell_filter, raw, origin=[0] * 4)

def test_basic():
    """Test Maxwell filter basic version."""
    # Load testing data (raw, SSS std origin, SSS non-standard origin)
    raw = read_crop(raw_fname, (0., 1.))
    raw_err = read_crop(raw_fname).apply_proj()
    raw_erm = read_crop(erm_fname)
    assert_raises(RuntimeError, maxwell_filter, raw_err)
    assert_raises(TypeError, maxwell_filter, 1.)  # not a raw
    assert_raises(ValueError, maxwell_filter, raw, int_order=20)  # too many

    n_int_bases = int_order ** 2 + 2 * int_order
    n_ext_bases = ext_order ** 2 + 2 * ext_order
    nbases = n_int_bases + n_ext_bases

    # Check number of bases computed correctly
    assert_equal(_get_n_moments([int_order, ext_order]).sum(), nbases)

    # Test SSS computation at the standard head origin
    assert_equal(len(raw.info['projs']), 12)  # 11 MEG projs + 1 AVG EEG
    raw_sss = maxwell_filter(raw, origin=mf_head_origin, regularize=None,
                             bad_condition='ignore')
    assert_equal(len(raw_sss.info['projs']), 1)  # avg EEG
    assert_equal(raw_sss.info['projs'][0]['desc'], 'Average EEG reference')
    assert_meg_snr(raw_sss, read_crop(sss_std_fname), 200., 1000.)
    py_cal = raw_sss.info['proc_history'][0]['max_info']['sss_cal']
    assert_equal(len(py_cal), 0)
    py_ctc = raw_sss.info['proc_history'][0]['max_info']['sss_ctc']
    assert_equal(len(py_ctc), 0)
    py_st = raw_sss.info['proc_history'][0]['max_info']['max_st']
    assert_equal(len(py_st), 0)
    assert_raises(RuntimeError, maxwell_filter, raw_sss)

    # Test SSS computation at non-standard head origin
    raw_sss = maxwell_filter(raw, origin=[0., 0.02, 0.02], regularize=None,
                             bad_condition='ignore')
    assert_meg_snr(raw_sss, read_crop(sss_nonstd_fname), 250., 700.)

    # Test SSS computation at device origin
    sss_erm_std = read_crop(sss_erm_std_fname)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg',
                             origin=mf_meg_origin, regularize=None,
                             bad_condition='ignore')
    assert_meg_snr(raw_sss, sss_erm_std, 100., 900.)
    for key in ('job', 'frame'):
        vals = [x.info['proc_history'][0]['max_info']['sss_info'][key]
                for x in [raw_sss, sss_erm_std]]
        assert_equal(vals[0], vals[1])

    # Check against SSS functions from proc_history
    sss_info = raw_sss.info['proc_history'][0]['max_info']
    assert_equal(_get_n_moments(int_order),
                 proc_history._get_sss_rank(sss_info))

    # Degenerate cases
    assert_raises(ValueError, maxwell_filter, raw, coord_frame='foo')
    assert_raises(ValueError, maxwell_filter, raw, origin='foo')
    assert_raises(ValueError, maxwell_filter, raw, origin=[0] * 4)
    assert_raises(ValueError, maxwell_filter, raw, mag_scale='foo')

def test_rank():
    """Test cov rank estimation."""
    # Test that our rank estimation works properly on a simple case
    evoked = read_evokeds(ave_fname, condition=0, baseline=(None, 0),
                          proj=False)
    cov = read_cov(cov_fname)
    ch_names = [ch for ch in evoked.info['ch_names'] if '053' not in ch and
                ch.startswith('EEG')]
    cov = prepare_noise_cov(cov, evoked.info, ch_names, None)
    assert_equal(cov['eig'][0], 0.)  # avg projector should set this to zero
    assert_true((cov['eig'][1:] > 0).all())  # all else should be > 0

    # Now do some more comprehensive tests
    raw_sample = read_raw_fif(raw_fname)

    raw_sss = read_raw_fif(hp_fif_fname)
    raw_sss.add_proj(compute_proj_raw(raw_sss))

    cov_sample = compute_raw_covariance(raw_sample)
    cov_sample_proj = compute_raw_covariance(
        raw_sample.copy().apply_proj())

    cov_sss = compute_raw_covariance(raw_sss)
    cov_sss_proj = compute_raw_covariance(
        raw_sss.copy().apply_proj())

    picks_all_sample = pick_types(raw_sample.info, meg=True, eeg=True)
    picks_all_sss = pick_types(raw_sss.info, meg=True, eeg=True)

    info_sample = pick_info(raw_sample.info, picks_all_sample)
    picks_stack_sample = [('eeg', pick_types(info_sample, meg=False,
                                             eeg=True))]
    picks_stack_sample += [('meg', pick_types(info_sample, meg=True))]
    picks_stack_sample += [('all',
                            pick_types(info_sample, meg=True, eeg=True))]

    info_sss = pick_info(raw_sss.info, picks_all_sss)
    picks_stack_somato = [('eeg', pick_types(info_sss, meg=False, eeg=True))]
    picks_stack_somato += [('meg', pick_types(info_sss, meg=True))]
    picks_stack_somato += [('all',
                            pick_types(info_sss, meg=True, eeg=True))]

    iter_tests = list(itt.product(
        [(cov_sample, picks_stack_sample, info_sample),
         (cov_sample_proj, picks_stack_sample, info_sample),
         (cov_sss, picks_stack_somato, info_sss),
         (cov_sss_proj, picks_stack_somato, info_sss)],  # sss
        [dict(mag=1e15, grad=1e13, eeg=1e6)]
    ))

    for (cov, picks_list, this_info), scalings in iter_tests:
        for ch_type, picks in picks_list:

            this_very_info = pick_info(this_info, picks)

            # compute subset of projs
            this_projs = [c['active'] and
                          len(set(c['data']['col_names'])
                              .intersection(set(this_very_info['ch_names']))) >
                          0 for c in cov['projs']]
            n_projs = sum(this_projs)

            # count channel types
            ch_types = [channel_type(this_very_info, idx)
                        for idx in range(len(picks))]
            n_eeg, n_mag, n_grad = [ch_types.count(k) for k in
                                    ['eeg', 'mag', 'grad']]
            n_meg = n_mag + n_grad
            if ch_type in ('all', 'eeg'):
                n_projs_eeg = 1
            else:
                n_projs_eeg = 0

            # check sss
            if len(this_very_info['proc_history']) > 0:
                mf = this_very_info['proc_history'][0]['max_info']
                n_free = _get_sss_rank(mf)
                if 'mag' not in ch_types and 'grad' not in ch_types:
                    n_free = 0
                # - n_projs XXX clarify
                expected_rank = n_free + n_eeg
                if n_projs > 0 and ch_type in ('all', 'eeg'):
                    expected_rank -= n_projs_eeg
            else:
                expected_rank = n_meg + n_eeg - n_projs

            C = cov['data'][np.ix_(picks, picks)]
            est_rank = _estimate_rank_meeg_cov(C, this_very_info,
                                               scalings=scalings)

            assert_equal(expected_rank, est_rank)

def test_basic():
    """Test Maxwell filter basic version."""
    # Load testing data (raw, SSS std origin, SSS non-standard origin)
    raw = read_crop(raw_fname, (0., 1.))
    raw_err = read_crop(raw_fname).apply_proj()
    raw_erm = read_crop(erm_fname)
    pytest.raises(RuntimeError, maxwell_filter, raw_err)
    pytest.raises(TypeError, maxwell_filter, 1.)  # not a raw
    pytest.raises(ValueError, maxwell_filter, raw, int_order=20)  # too many

    n_int_bases = int_order ** 2 + 2 * int_order
    n_ext_bases = ext_order ** 2 + 2 * ext_order
    nbases = n_int_bases + n_ext_bases

    # Check number of bases computed correctly
    assert_equal(_get_n_moments([int_order, ext_order]).sum(), nbases)

    # Test SSS computation at the standard head origin
    assert_equal(len(raw.info['projs']), 12)  # 11 MEG projs + 1 AVG EEG
    with use_coil_def(elekta_def_fname):
        raw_sss = maxwell_filter(raw, origin=mf_head_origin, regularize=None,
                                 bad_condition='ignore')
    assert_equal(len(raw_sss.info['projs']), 1)  # avg EEG
    assert_equal(raw_sss.info['projs'][0]['desc'], 'Average EEG reference')
    assert_meg_snr(raw_sss, read_crop(sss_std_fname), 200., 1000.)
    py_cal = raw_sss.info['proc_history'][0]['max_info']['sss_cal']
    assert_equal(len(py_cal), 0)
    py_ctc = raw_sss.info['proc_history'][0]['max_info']['sss_ctc']
    assert_equal(len(py_ctc), 0)
    py_st = raw_sss.info['proc_history'][0]['max_info']['max_st']
    assert_equal(len(py_st), 0)
    pytest.raises(RuntimeError, maxwell_filter, raw_sss)

    # Test SSS computation at non-standard head origin
    with use_coil_def(elekta_def_fname):
        raw_sss = maxwell_filter(raw, origin=[0., 0.02, 0.02], regularize=None,
                                 bad_condition='ignore')
    assert_meg_snr(raw_sss, read_crop(sss_nonstd_fname), 250., 700.)

    # Test SSS computation at device origin
    sss_erm_std = read_crop(sss_erm_std_fname)
    raw_sss = maxwell_filter(raw_erm, coord_frame='meg',
                             origin=mf_meg_origin, regularize=None,
                             bad_condition='ignore')
    assert_meg_snr(raw_sss, sss_erm_std, 70., 260.)
    for key in ('job', 'frame'):
        vals = [x.info['proc_history'][0]['max_info']['sss_info'][key]
                for x in [raw_sss, sss_erm_std]]
        assert_equal(vals[0], vals[1])

    # Two equivalent things: at device origin in device coords (0., 0., 0.)
    # and at device origin at head coords info['dev_head_t'][:3, 3]
    raw_sss_meg = maxwell_filter(
        raw, coord_frame='meg', origin=(0., 0., 0.))
    raw_sss_head = maxwell_filter(
        raw, origin=raw.info['dev_head_t']['trans'][:3, 3])
    assert_meg_snr(raw_sss_meg, raw_sss_head, 100., 900.)

    # Check against SSS functions from proc_history
    sss_info = raw_sss.info['proc_history'][0]['max_info']
    assert_equal(_get_n_moments(int_order),
                 proc_history._get_sss_rank(sss_info))

    # Degenerate cases
    pytest.raises(ValueError, maxwell_filter, raw, coord_frame='foo')
    pytest.raises(ValueError, maxwell_filter, raw, origin='foo')
    pytest.raises(ValueError, maxwell_filter, raw, origin=[0] * 4)
    pytest.raises(ValueError, maxwell_filter, raw, mag_scale='foo')
    raw_missing = raw.copy().load_data()
    raw_missing.info['bads'] = ['MEG0111']
    raw_missing.pick_types(meg=True)  # will be missing the bad
    maxwell_filter(raw_missing)
    with pytest.warns(RuntimeWarning, match='not in data'):
        maxwell_filter(raw_missing, calibration=fine_cal_fname)

def test_rank():
    """Test cov rank estimation"""
    raw_sample = Raw(raw_fname)

    raw_sss = Raw(hp_fif_fname)
    raw_sss.add_proj(compute_proj_raw(raw_sss))

    cov_sample = compute_raw_data_covariance(raw_sample)
    cov_sample_proj = compute_raw_data_covariance(
        raw_sample.copy().apply_proj())

    cov_sss = compute_raw_data_covariance(raw_sss)
    cov_sss_proj = compute_raw_data_covariance(
        raw_sss.copy().apply_proj())

    picks_all_sample = pick_types(raw_sample.info, meg=True, eeg=True)
    picks_all_sss = pick_types(raw_sss.info, meg=True, eeg=True)

    info_sample = pick_info(raw_sample.info, picks_all_sample)
    picks_stack_sample = [('eeg', pick_types(info_sample, meg=False,
                                             eeg=True))]
    picks_stack_sample += [('meg', pick_types(info_sample, meg=True))]
    picks_stack_sample += [('all',
                            pick_types(info_sample, meg=True, eeg=True))]

    info_sss = pick_info(raw_sss.info, picks_all_sss)
    picks_stack_somato = [('eeg', pick_types(info_sss, meg=False, eeg=True))]
    picks_stack_somato += [('meg', pick_types(info_sss, meg=True))]
    picks_stack_somato += [('all',
                            pick_types(info_sss, meg=True, eeg=True))]

    iter_tests = list(itt.product(
        [(cov_sample, picks_stack_sample, info_sample),
         (cov_sample_proj, picks_stack_sample, info_sample),
         (cov_sss, picks_stack_somato, info_sss),
         (cov_sss_proj, picks_stack_somato, info_sss)],  # sss
        [dict(mag=1e15, grad=1e13, eeg=1e6)]
    ))

    for (cov, picks_list, this_info), scalings in iter_tests:
        for ch_type, picks in picks_list:

            this_very_info = pick_info(this_info, picks)

            # compute subset of projs
            this_projs = [c['active'] and
                          len(set(c['data']['col_names'])
                              .intersection(set(this_very_info['ch_names']))) >
                          0 for c in cov['projs']]
            n_projs = sum(this_projs)

            # count channel types
            ch_types = [channel_type(this_very_info, idx)
                        for idx in range(len(picks))]
            n_eeg, n_mag, n_grad = [ch_types.count(k) for k in
                                    ['eeg', 'mag', 'grad']]
            n_meg = n_mag + n_grad
            if ch_type in ('all', 'eeg'):
                n_projs_eeg = 1
            else:
                n_projs_eeg = 0

            # check sss
            if 'proc_history' in this_very_info:
                mf = this_very_info['proc_history'][0]['max_info']
                n_free = _get_sss_rank(mf)
                if 'mag' not in ch_types and 'grad' not in ch_types:
                    n_free = 0
                # - n_projs XXX clarify
                expected_rank = n_free + n_eeg
                if n_projs > 0 and ch_type in ('all', 'eeg'):
                    expected_rank -= n_projs_eeg
            else:
                expected_rank = n_meg + n_eeg - n_projs

            C = cov['data'][np.ix_(picks, picks)]
            est_rank = _estimate_rank_meeg_cov(C, this_very_info,
                                               scalings=scalings)

            assert_equal(expected_rank, est_rank)

def test_maxwell_filter():
    """Test multipolar moment and Maxwell filter"""

    # TODO: Future tests integrate with mne/io/tests/test_proc_history

    # Load testing data (raw, SSS std origin, SSS non-standard origin)
    with warnings.catch_warnings(record=True):  # maxshield
        raw = Raw(raw_fname, preload=False, proj=False,
                  allow_maxshield=True).crop(0., 1., False)
    raw.load_data()
    with warnings.catch_warnings(record=True):  # maxshield, naming
        sss_std = Raw(sss_std_fname, preload=True, proj=False,
                      allow_maxshield=True)
        sss_nonStd = Raw(sss_nonstd_fname, preload=True, proj=False,
                         allow_maxshield=True)
        raw_err = Raw(raw_fname, preload=False, proj=True,
                      allow_maxshield=True).crop(0., 0.1, False)
    assert_raises(RuntimeError, maxwell.maxwell_filter, raw_err)

    # Create coils
    all_coils, _, _, meg_info = _prep_meg_channels(raw.info, ignore_ref=True,
                                                   elekta_defs=True)
    picks = [raw.info['ch_names'].index(ch) for ch in [coil['chname']
                                                       for coil in all_coils]]
    coils = [all_coils[ci] for ci in picks]
    ncoils = len(coils)

    int_order, ext_order = 8, 3
    n_int_bases = int_order ** 2 + 2 * int_order
    n_ext_bases = ext_order ** 2 + 2 * ext_order
    nbases = n_int_bases + n_ext_bases

    # Check number of bases computed correctly
    assert_equal(maxwell.get_num_moments(int_order, ext_order), nbases)

    # Check multipolar moment basis set
    S_in, S_out = maxwell._sss_basis(origin=np.array([0, 0, 40]), coils=coils,
                                     int_order=int_order, ext_order=ext_order)
    assert_equal(S_in.shape, (ncoils, n_int_bases), 'S_in has incorrect shape')
    assert_equal(S_out.shape, (ncoils, n_ext_bases),
                 'S_out has incorrect shape')

    # Test sss computation at the standard head origin
    raw_sss = maxwell.maxwell_filter(raw, origin=[0., 0., 40.],
                                     int_order=int_order, ext_order=ext_order)

    assert_array_almost_equal(raw_sss._data[picks, :], sss_std._data[picks, :],
                              decimal=11, err_msg='Maxwell filtered data at '
                              'standard origin incorrect.')

    # Confirm SNR is above 100
    bench_rms = np.sqrt(np.mean(sss_std._data[picks, :] ** 2, axis=1))
    error = raw_sss._data[picks, :] - sss_std._data[picks, :]
    error_rms = np.sqrt(np.mean(error ** 2, axis=1))
    assert_true(np.mean(bench_rms / error_rms) > 1000, 'SNR < 1000')

    # Test sss computation at non-standard head origin
    raw_sss = maxwell.maxwell_filter(raw, origin=[0., 20., 20.],
                                     int_order=int_order, ext_order=ext_order)
    assert_array_almost_equal(raw_sss._data[picks, :],
                              sss_nonStd._data[picks, :], decimal=11,
                              err_msg='Maxwell filtered data at non-std '
                              'origin incorrect.')
    # Confirm SNR is above 100
    bench_rms = np.sqrt(np.mean(sss_nonStd._data[picks, :] ** 2, axis=1))
    error = raw_sss._data[picks, :] - sss_nonStd._data[picks, :]
    error_rms = np.sqrt(np.mean(error ** 2, axis=1))
    assert_true(np.mean(bench_rms / error_rms) > 1000, 'SNR < 1000')

    # Check against SSS functions from proc_history
    sss_info = raw_sss.info['proc_history'][0]['max_info']
    assert_equal(maxwell.get_num_moments(int_order, 0),
                 proc_history._get_sss_rank(sss_info))

def test_maxfilter_get_rank():
    """Test maxfilter rank lookup."""
    mf = read_info(raw_fname)['proc_history'][0]['max_info']
    rank1 = mf['sss_info']['nfree']
    rank2 = _get_sss_rank(mf)
    assert_equal(rank1, rank2)