Python ICA.pick_sources_raw方法代码示例

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import pick_sources_raw [as 别名]
def test_ica_twice():
    """Test running ICA twice"""
    raw = io.Raw(raw_fname, preload=True).crop(0, stop, False).crop(1.5)
    picks = pick_types(raw.info, meg='grad', exclude='bads')
    n_components = 0.9
    max_pca_components = None
    n_pca_components = 1.1
    with warnings.catch_warnings(record=True):
        ica1 = ICA(n_components=n_components, max_pca_components=max_pca_components,
                   n_pca_components=n_pca_components, random_state=0)

        ica1.decompose_raw(raw, picks=picks, decim=3)
        raw_new = ica1.pick_sources_raw(raw, n_pca_components=n_pca_components)
        ica2 = ICA(n_components=n_components, max_pca_components=max_pca_components,
                   n_pca_components=1.0, random_state=0)
        ica2.decompose_raw(raw_new, picks=picks, decim=3)
        assert_equal(ica1.n_components_, ica2.n_components_)

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import pick_sources_raw [as 别名]
def test_ica_full_data_recovery():
    """Test recovery of full data when no source is rejected"""
    # Most basic recovery
    raw = io.Raw(raw_fname, preload=True).crop(0, stop, False).crop(1.5)
    events = read_events(event_name)
    picks = pick_types(raw.info, meg=True, stim=False, ecg=False,
                       eog=False, exclude='bads')
    epochs = Epochs(raw, events[:4], event_id, tmin, tmax, picks=picks,
                    baseline=(None, 0), preload=True)
    n_channels = 5
    data = raw._data[:n_channels].copy()
    data_epochs = epochs.get_data()
    for n_components, n_pca_components, ok in [(2, n_channels, True),
                                               (2, n_channels // 2, False)]:
        ica = ICA(n_components=n_components,
                  max_pca_components=n_pca_components,
                  n_pca_components=n_pca_components)
        ica.decompose_raw(raw, picks=list(range(n_channels)))
        raw2 = ica.pick_sources_raw(raw, exclude=[])
        if ok:
            assert_allclose(data[:n_channels], raw2._data[:n_channels],
                            rtol=1e-10, atol=1e-15)
        else:
            diff = np.abs(data[:n_channels] - raw2._data[:n_channels])
            assert_true(np.max(diff) > 1e-14)

        ica = ICA(n_components=n_components,
                  max_pca_components=n_pca_components,
                  n_pca_components=n_pca_components)
        ica.decompose_epochs(epochs, picks=list(range(n_channels)))
        epochs2 = ica.pick_sources_epochs(epochs, exclude=[])
        data2 = epochs2.get_data()[:, :n_channels]
        if ok:
            assert_allclose(data_epochs[:, :n_channels], data2,
                            rtol=1e-10, atol=1e-15)
        else:
            diff = np.abs(data_epochs[:, :n_channels] - data2)
            assert_true(np.max(diff) > 1e-14)

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import pick_sources_raw [as 别名]
def test_ica_additional():
    """Test additional ICA functionality
    """
    stop2 = 500
    raw = io.Raw(raw_fname, preload=True).crop(0, stop, False).crop(1.5)
    picks = pick_types(raw.info, meg=True, stim=False, ecg=False,
                       eog=False, exclude='bads')
    test_cov = read_cov(test_cov_name)
    events = read_events(event_name)
    picks = pick_types(raw.info, meg=True, stim=False, ecg=False,
                       eog=False, exclude='bads')
    epochs = Epochs(raw, events[:4], event_id, tmin, tmax, picks=picks,
                    baseline=(None, 0), preload=True)
    # for testing eog functionality
    picks2 = pick_types(raw.info, meg=True, stim=False, ecg=False,
                        eog=True, exclude='bads')
    epochs_eog = Epochs(raw, events[:4], event_id, tmin, tmax, picks=picks2,
                        baseline=(None, 0), preload=True)

    test_cov2 = deepcopy(test_cov)
    ica = ICA(noise_cov=test_cov2, n_components=3, max_pca_components=4,
              n_pca_components=4)
    assert_true(ica.info is None)
    ica.decompose_raw(raw, picks[:5])
    assert_true(isinstance(ica.info, Info))
    assert_true(ica.n_components_ < 5)

    ica = ICA(n_components=3, max_pca_components=4,
              n_pca_components=4)
    assert_raises(RuntimeError, ica.save, '')
    ica.decompose_raw(raw, picks=None, start=start, stop=stop2)

    # test warnings on bad filenames
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter('always')
        ica_badname = op.join(op.dirname(tempdir), 'test-bad-name.fif.gz')
        ica.save(ica_badname)
        read_ica(ica_badname)
    assert_true(len(w) == 2)

    # test decim
    ica = ICA(n_components=3, max_pca_components=4,
              n_pca_components=4)
    raw_ = raw.copy()
    for _ in range(3):
        raw_.append(raw_)
    n_samples = raw_._data.shape[1]
    ica.decompose_raw(raw, picks=None, decim=3)
    assert_true(raw_._data.shape[1], n_samples)

    # test expl var
    ica = ICA(n_components=1.0, max_pca_components=4,
              n_pca_components=4)
    ica.decompose_raw(raw, picks=None, decim=3)
    assert_true(ica.n_components_ == 4)

    # epochs extraction from raw fit
    assert_raises(RuntimeError, ica.get_sources_epochs, epochs)
    # test reading and writing
    test_ica_fname = op.join(op.dirname(tempdir), 'test-ica.fif')
    for cov in (None, test_cov):
        ica = ICA(noise_cov=cov, n_components=2, max_pca_components=4,
                  n_pca_components=4)
        with warnings.catch_warnings(record=True):  # ICA does not converge
            ica.decompose_raw(raw, picks=picks, start=start, stop=stop2)
        sources = ica.get_sources_epochs(epochs)
        assert_true(ica.mixing_matrix_.shape == (2, 2))
        assert_true(ica.unmixing_matrix_.shape == (2, 2))
        assert_true(ica.pca_components_.shape == (4, len(picks)))
        assert_true(sources.shape[1] == ica.n_components_)

        for exclude in [[], [0]]:
            ica.exclude = [0]
            ica.save(test_ica_fname)
            ica_read = read_ica(test_ica_fname)
            assert_true(ica.exclude == ica_read.exclude)
            # test pick merge -- add components
            ica.pick_sources_raw(raw, exclude=[1])
            assert_true(ica.exclude == [0, 1])
            #                 -- only as arg
            ica.exclude = []
            ica.pick_sources_raw(raw, exclude=[0, 1])
            assert_true(ica.exclude == [0, 1])
            #                 -- remove duplicates
            ica.exclude += [1]
            ica.pick_sources_raw(raw, exclude=[0, 1])
            assert_true(ica.exclude == [0, 1])

            # test basic include
            ica.exclude = []
            ica.pick_sources_raw(raw, include=[1])

            ica_raw = ica.sources_as_raw(raw)
            assert_true(ica.exclude == [ica_raw.ch_names.index(e) for e in
                                        ica_raw.info['bads']])

        # test filtering
        d1 = ica_raw._data[0].copy()
        with warnings.catch_warnings(record=True):  # dB warning
            ica_raw.filter(4, 20)
#.........这里部分代码省略.........

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import pick_sources_raw [as 别名]
    
    # Select ICs which can be attributed to cardiac artifacts
    Brain_sources = pk_max > 0.1  # bool array, get the prominient components related with trigger
    Brain_ind = np.where(Brain_sources)[0].tolist() # indices
    #skip the null idx related with response
    Brain_idx += (Brain_ind)#Get the obvious sources related
    if Brain_ind == []:
    #plot_phase = False
    #if plot_phase == False: 
      continue
  
Brain_idx = list(set(Brain_idx))
print '%s has been identified as Brain response' %(Brain_idx)

#Recompose raw_cle data include the response
raw_brain = ica.pick_sources_raw(raw, include=Brain_idx, n_pca_components=n_pca_components)
raw_brain.save(subject_path+'/MEG/raw_%s.fif' %(raw_basename), overwrite=True)#Get the raw data related with brain response

#========Compare the brain response between raw before and after ctps=====
if trigger == 'resp':
    tmin, tmax = -0.3, 0.3
    events = mne.find_events(raw, stim_channel=res_ch_name)
elif trigger == 'stim':
    tmin, tmax = -0.2, 0.4
    events = mne.find_events(raw, stim_channel=sti_ch_name)
picks_bef = mne.fiff.pick_types(raw.info, meg=True, exclude='bads')
picks_aft = mne.fiff.pick_types(raw_brain.info, meg=True, exclude='bads')
evoked_before_ctps = mne.Epochs(raw, events, 1, tmin, tmax, proj=False, picks=picks_bef,   preload=False, reject=None).average()
evoked_after_ctps = mne.Epochs(raw_brain, events, 1, tmin, tmax, proj=False, picks=picks_aft,   preload=False, reject=None).average()
evoked_after_ctps.save(subject_path+'/MEG/ave_%s.fif' %(raw_basename))
times = evoked_before_ctps.times

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import pick_sources_raw [as 别名]
def test_ica_core():
    """Test ICA on raw and epochs
    """
    # setup parameter
    # XXX. The None cases helped revealing bugs but are time consuming.
    noise_cov = [None, test_cov]
    # removed None cases to speed up...
    n_components = [3, 1.0]  # for future dbg add cases
    max_pca_components = [4]
    picks_ = [picks]
    iter_ica_params = product(noise_cov, n_components, max_pca_components,
                              picks_)

    # # test init catchers
    assert_raises(ValueError, ICA, n_components=3, max_pca_components=2)
    assert_raises(ValueError, ICA, n_components=1.3, max_pca_components=2)

    # test essential core functionality
    for n_cov, n_comp, max_n, pcks in iter_ica_params:
      # Test ICA raw
        ica = ICA(noise_cov=n_cov, n_components=n_comp,
                  max_pca_components=max_n, n_pca_components=max_n,
                  random_state=0)

        print ica  # to test repr

        # test fit checker
        assert_raises(RuntimeError, ica.get_sources_raw, raw)
        assert_raises(RuntimeError, ica.get_sources_epochs, epochs)

        # test decomposition
        ica.decompose_raw(raw, picks=pcks, start=start, stop=stop)
        print ica  # to test repr
        # test re-init exception
        assert_raises(RuntimeError, ica.decompose_raw, raw, picks=picks)

        sources = ica.get_sources_raw(raw)
        assert_true(sources.shape[0] == ica.n_components_)

        # test preload filter
        raw3 = raw.copy()
        raw3._preloaded = False
        assert_raises(ValueError, ica.pick_sources_raw, raw3,
                      include=[1, 2])

        for excl, incl in (([], []), ([], [1, 2]), ([1, 2], [])):
            raw2 = ica.pick_sources_raw(raw, exclude=excl, include=incl,
                                        copy=True)

            assert_array_almost_equal(raw2[:, :][1], raw[:, :][1])

        #######################################################################
        # test epochs decomposition

        # test re-init exception
        assert_raises(RuntimeError, ica.decompose_epochs, epochs, picks=picks)
        ica = ICA(noise_cov=n_cov, n_components=n_comp,
                  max_pca_components=max_n, n_pca_components=max_n,
                  random_state=0)

        ica.decompose_epochs(epochs, picks=picks)
        print ica  # to test repr
        # test pick block after epochs fit
        assert_raises(ValueError, ica.pick_sources_raw, raw)

        sources = ica.get_sources_epochs(epochs)
        assert_true(sources.shape[1] == ica.n_components_)

        assert_raises(ValueError, ica.find_sources_epochs, epochs,
                      target=np.arange(1))

        # test preload filter
        epochs3 = epochs.copy()
        epochs3.preload = False
        assert_raises(ValueError, ica.pick_sources_epochs, epochs3,
                      include=[1, 2])

        # test source picking
        for excl, incl in (([], []), ([], [1, 2]), ([1, 2], [])):
            epochs2 = ica.pick_sources_epochs(epochs, exclude=excl,
                                      include=incl, copy=True)

            assert_array_almost_equal(epochs2.get_data(),
                                      epochs.get_data())

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import pick_sources_raw [as 别名]
def test_ica_additional():
    """Test additional functionality
    """
    stop2 = 500

    test_cov2 = deepcopy(test_cov)
    ica = ICA(noise_cov=test_cov2, n_components=3, max_pca_components=4,
              n_pca_components=4)
    ica.decompose_raw(raw, picks[:5])
    assert_true(ica.n_components_ < 5)

    ica = ICA(n_components=3, max_pca_components=4,
              n_pca_components=4)
    assert_raises(RuntimeError, ica.save, '')
    ica.decompose_raw(raw, picks=None, start=start, stop=stop2)

    # epochs extraction from raw fit
    assert_raises(RuntimeError, ica.get_sources_epochs, epochs)

    # test reading and writing
    test_ica_fname = op.join(op.dirname(tempdir), 'ica_test.fif')
    for cov in (None, test_cov):
        ica = ICA(noise_cov=cov, n_components=3, max_pca_components=4,
                  n_pca_components=4)
        ica.decompose_raw(raw, picks=picks, start=start, stop=stop2)
        sources = ica.get_sources_epochs(epochs)
        assert_true(sources.shape[1] == ica.n_components_)

        for exclude in [[], [0]]:
            ica.exclude = [0]
            ica.save(test_ica_fname)
            ica_read = read_ica(test_ica_fname)
            assert_true(ica.exclude == ica_read.exclude)
            # test pick merge -- add components
            ica.pick_sources_raw(raw, exclude=[1])
            assert_true(ica.exclude == [0, 1])
            #                 -- only as arg
            ica.exclude = []
            ica.pick_sources_raw(raw, exclude=[0, 1])
            assert_true(ica.exclude == [0, 1])
            #                 -- remove duplicates
            ica.exclude += [1]
            ica.pick_sources_raw(raw, exclude=[0, 1])
            assert_true(ica.exclude == [0, 1])

            ica_raw = ica.sources_as_raw(raw)
            assert_true(ica.exclude == [ica.ch_names.index(e) for e in
                                        ica_raw.info['bads']])

        ica.n_pca_components = 2
        ica.save(test_ica_fname)
        ica_read = read_ica(test_ica_fname)
        assert_true(ica.n_pca_components ==
                    ica_read.n_pca_components)
        ica.n_pca_components = 4
        ica_read.n_pca_components = 4

        ica.exclude = []
        ica.save(test_ica_fname)
        ica_read = read_ica(test_ica_fname)

        assert_true(ica.ch_names == ica_read.ch_names)

        assert_true(np.allclose(ica.mixing_matrix_, ica_read.mixing_matrix_,
                                rtol=1e-16, atol=1e-32))
        assert_array_equal(ica.pca_components_,
                           ica_read.pca_components_)
        assert_array_equal(ica.pca_mean_, ica_read.pca_mean_)
        assert_array_equal(ica.pca_explained_variance_,
                           ica_read.pca_explained_variance_)
        assert_array_equal(ica._pre_whitener, ica_read._pre_whitener)

        # assert_raises(RuntimeError, ica_read.decompose_raw, raw)
        sources = ica.get_sources_raw(raw)
        sources2 = ica_read.get_sources_raw(raw)
        assert_array_almost_equal(sources, sources2)

        _raw1 = ica.pick_sources_raw(raw, exclude=[1])
        _raw2 = ica_read.pick_sources_raw(raw, exclude=[1])
        assert_array_almost_equal(_raw1[:, :][0], _raw2[:, :][0])

    os.remove(test_ica_fname)
    # score funcs raw, with catch since "ties preclude exact" warning
    # XXX this should be fixed by a future PR...
    with warnings.catch_warnings(True) as w:
        sfunc_test = [ica.find_sources_raw(raw, target='EOG 061',
                score_func=n, start=0, stop=10)
                for n, f in score_funcs.items()]
    # score funcs raw

    # check lenght of scores
    [assert_true(ica.n_components_ == len(scores)) for scores in sfunc_test]

    # check univariate stats
    scores = ica.find_sources_raw(raw, score_func=stats.skew)
    # check exception handling
    assert_raises(ValueError, ica.find_sources_raw, raw,
                  target=np.arange(1))

    ## score funcs epochs ##
#.........这里部分代码省略.........

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import pick_sources_raw [as 别名]
def test_ica_additional():
    """Test additional functionality
    """
    stop2 = 500

    test_cov2 = deepcopy(test_cov)
    ica = ICA(noise_cov=test_cov2, n_components=3, max_pca_components=4,
              n_pca_components=4)
    assert_true(ica.info is None)
    ica.decompose_raw(raw, picks[:5])
    assert_true(isinstance(ica.info, Info))
    assert_true(ica.n_components_ < 5)

    ica = ICA(n_components=3, max_pca_components=4,
              n_pca_components=4)
    assert_raises(RuntimeError, ica.save, '')
    ica.decompose_raw(raw, picks=None, start=start, stop=stop2)

    # epochs extraction from raw fit
    assert_raises(RuntimeError, ica.get_sources_epochs, epochs)

    # test reading and writing
    test_ica_fname = op.join(op.dirname(tempdir), 'ica_test.fif')
    for cov in (None, test_cov):
        ica = ICA(noise_cov=cov, n_components=3, max_pca_components=4,
                  n_pca_components=4)
        ica.decompose_raw(raw, picks=picks, start=start, stop=stop2)
        sources = ica.get_sources_epochs(epochs)
        assert_true(sources.shape[1] == ica.n_components_)

        for exclude in [[], [0]]:
            ica.exclude = [0]
            ica.save(test_ica_fname)
            ica_read = read_ica(test_ica_fname)
            assert_true(ica.exclude == ica_read.exclude)
            # test pick merge -- add components
            ica.pick_sources_raw(raw, exclude=[1])
            assert_true(ica.exclude == [0, 1])
            #                 -- only as arg
            ica.exclude = []
            ica.pick_sources_raw(raw, exclude=[0, 1])
            assert_true(ica.exclude == [0, 1])
            #                 -- remove duplicates
            ica.exclude += [1]
            ica.pick_sources_raw(raw, exclude=[0, 1])
            assert_true(ica.exclude == [0, 1])

            ica_raw = ica.sources_as_raw(raw)
            assert_true(ica.exclude == [ica_raw.ch_names.index(e) for e in
                                        ica_raw.info['bads']])

        ica.n_pca_components = 2
        ica.save(test_ica_fname)
        ica_read = read_ica(test_ica_fname)
        assert_true(ica.n_pca_components ==
                    ica_read.n_pca_components)
        ica.n_pca_components = 4
        ica_read.n_pca_components = 4

        ica.exclude = []
        ica.save(test_ica_fname)
        ica_read = read_ica(test_ica_fname)

        assert_true(ica.ch_names == ica_read.ch_names)
        assert_true(isinstance(ica_read.info, Info))  # XXX improve later
        assert_true(np.allclose(ica.mixing_matrix_, ica_read.mixing_matrix_,
                                rtol=1e-16, atol=1e-32))
        assert_array_equal(ica.pca_components_,
                           ica_read.pca_components_)
        assert_array_equal(ica.pca_mean_, ica_read.pca_mean_)
        assert_array_equal(ica.pca_explained_variance_,
                           ica_read.pca_explained_variance_)
        assert_array_equal(ica._pre_whitener, ica_read._pre_whitener)

        # assert_raises(RuntimeError, ica_read.decompose_raw, raw)
        sources = ica.get_sources_raw(raw)
        sources2 = ica_read.get_sources_raw(raw)
        assert_array_almost_equal(sources, sources2)

        _raw1 = ica.pick_sources_raw(raw, exclude=[1])
        _raw2 = ica_read.pick_sources_raw(raw, exclude=[1])
        assert_array_almost_equal(_raw1[:, :][0], _raw2[:, :][0])

    os.remove(test_ica_fname)
    # check scrore funcs
    for name, func in score_funcs.items():
        if name in score_funcs_unsuited:
            continue
        scores = ica.find_sources_raw(raw, target='EOG 061', score_func=func,
                                      start=0, stop=10)
        assert_true(ica.n_components_ == len(scores))

    # check univariate stats
    scores = ica.find_sources_raw(raw, score_func=stats.skew)
    # check exception handling
    assert_raises(ValueError, ica.find_sources_raw, raw,
                  target=np.arange(1))

    params = []
    params += [(None, -1, slice(2), [0, 1])]  # varicance, kurtosis idx params
#.........这里部分代码省略.........

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import pick_sources_raw [as 别名]
sources_ecg.filter(l_freq=8, h_freq=16, method='iir', n_jobs=4)

# Epochs at R peak onset, from ecg_eve.
ica_epochs_ecg = mne.Epochs(sources_ecg, ecg_eve, event_id=999, tmin=-0.5, tmax=0.5,
                        picks=picks, preload=True, proj=False)

# Compute phase values and statistics (significance values pK)
#phase_trial_ecg, pk_dyn_ecg, _ = compute_ctps(ica_epochs_ecg.get_data())
_ , pk_dyn_ecg, _ = compute_ctps(ica_epochs_ecg.get_data())

# Get kuiper maxima
pk_max = pk_dyn_ecg.max(axis=1)

# Select ICs which can be attributed to cardiac artifacts
reject_sources = pk_max > 0.20  # bool array
reject_idx_ecg = np.where(reject_sources)[0].tolist()  # indices
ica.exclude += reject_idx_ecg
print '%s has been identified as ecg components and excluded' %(reject_idx_ecg)


raw_cle = ica.pick_sources_raw(raw, n_pca_components=n_pca_components)
#if view_plots:
 #   import compare_raw_files
  #  compare_raw_files.compare_raw_files(raw, raw_cle)

#ica.save('ctps_cor_resp_013.fif')
raw_cle_fil = subject_path + '/MEG/%s_cle.fif' %(raw_basename)
raw_cle.save(raw_cle_fil, overwrite=True)#Get the cleaned raw data
import compare_raw_files
compare_raw_files.compare_raw_files(raw, raw_cle)