Python ICA.find_bads_eog方法代码示例

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]
def test_ica_labels():
    """Test ICA labels."""
    # The CTF data are uniquely well suited to testing the ICA.find_bads_
    # methods
    raw = read_raw_ctf(ctf_fname, preload=True)
    # derive reference ICA components and append them to raw
    icarf = ICA(n_components=2, random_state=0, max_iter=2, allow_ref_meg=True)
    with pytest.warns(UserWarning, match='did not converge'):
        icarf.fit(raw.copy().pick_types(meg=False, ref_meg=True))
    icacomps = icarf.get_sources(raw)
    # rename components so they are auto-detected by find_bads_ref
    icacomps.rename_channels({c: 'REF_' + c for c in icacomps.ch_names})
    # and add them to raw
    raw.add_channels([icacomps])
    # set the appropriate EEG channels to EOG and ECG
    raw.set_channel_types({'EEG057': 'eog', 'EEG058': 'eog', 'EEG059': 'ecg'})
    ica = ICA(n_components=4, random_state=0, max_iter=2, method='fastica')
    with pytest.warns(UserWarning, match='did not converge'):
        ica.fit(raw)

    ica.find_bads_eog(raw, l_freq=None, h_freq=None)
    picks = list(pick_types(raw.info, meg=False, eog=True))
    for idx, ch in enumerate(picks):
        assert '{}/{}/{}'.format('eog', idx, raw.ch_names[ch]) in ica.labels_
    assert 'eog' in ica.labels_
    for key in ('ecg', 'ref_meg', 'ecg/ECG-MAG'):
        assert key not in ica.labels_

    ica.find_bads_ecg(raw, l_freq=None, h_freq=None, method='correlation')
    picks = list(pick_types(raw.info, meg=False, ecg=True))
    for idx, ch in enumerate(picks):
        assert '{}/{}/{}'.format('ecg', idx, raw.ch_names[ch]) in ica.labels_
    for key in ('ecg', 'eog'):
        assert key in ica.labels_
    for key in ('ref_meg', 'ecg/ECG-MAG'):
        assert key not in ica.labels_

    ica.find_bads_ref(raw, l_freq=None, h_freq=None)
    picks = pick_channels_regexp(raw.ch_names, 'REF_ICA*')
    for idx, ch in enumerate(picks):
        assert '{}/{}/{}'.format('ref_meg', idx,
                                 raw.ch_names[ch]) in ica.labels_
    for key in ('ecg', 'eog', 'ref_meg'):
        assert key in ica.labels_
    assert 'ecg/ECG-MAG' not in ica.labels_

    ica.find_bads_ecg(raw, l_freq=None, h_freq=None)
    for key in ('ecg', 'eog', 'ref_meg', 'ecg/ECG-MAG'):
        assert key in ica.labels_

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]

#.........这里部分代码省略.........
                  target=np.arange(1))

    params = []
    params += [(None, -1, slice(2), [0, 1])]  # variance, kurtosis params
    params += [(None, 'MEG 1531')]  # ECG / EOG channel params
    for idx, ch_name in product(*params):
        ica.detect_artifacts(raw, start_find=0, stop_find=50, ecg_ch=ch_name,
                             eog_ch=ch_name, skew_criterion=idx,
                             var_criterion=idx, kurt_criterion=idx)

    # Make sure detect_artifacts marks the right components.
    # For int criterion, the doc says "E.g. range(2) would return the two
    # sources with the highest score". Assert that's what it does.
    # Only test for skew, since it's always the same code.
    ica.exclude = []
    ica.detect_artifacts(raw, start_find=0, stop_find=50, ecg_ch=None,
                         eog_ch=None, skew_criterion=0,
                         var_criterion=None, kurt_criterion=None)
    assert np.abs(scores[ica.exclude]) == np.max(np.abs(scores))

    evoked = epochs.average()
    evoked_data = evoked.data.copy()
    raw_data = raw[:][0].copy()
    epochs_data = epochs.get_data().copy()

    with pytest.warns(RuntimeWarning, match='longer'):
        idx, scores = ica.find_bads_ecg(raw, method='ctps')
    assert_equal(len(scores), ica.n_components_)
    with pytest.warns(RuntimeWarning, match='longer'):
        idx, scores = ica.find_bads_ecg(raw, method='correlation')
    assert_equal(len(scores), ica.n_components_)

    with pytest.warns(RuntimeWarning, match='longer'):
        idx, scores = ica.find_bads_eog(raw)
    assert_equal(len(scores), ica.n_components_)

    idx, scores = ica.find_bads_ecg(epochs, method='ctps')

    assert_equal(len(scores), ica.n_components_)
    pytest.raises(ValueError, ica.find_bads_ecg, epochs.average(),
                  method='ctps')
    pytest.raises(ValueError, ica.find_bads_ecg, raw,
                  method='crazy-coupling')

    with pytest.warns(RuntimeWarning, match='longer'):
        idx, scores = ica.find_bads_eog(raw)
    assert_equal(len(scores), ica.n_components_)

    raw.info['chs'][raw.ch_names.index('EOG 061') - 1]['kind'] = 202
    with pytest.warns(RuntimeWarning, match='longer'):
        idx, scores = ica.find_bads_eog(raw)
    assert (isinstance(scores, list))
    assert_equal(len(scores[0]), ica.n_components_)

    idx, scores = ica.find_bads_eog(evoked, ch_name='MEG 1441')
    assert_equal(len(scores), ica.n_components_)

    idx, scores = ica.find_bads_ecg(evoked, method='correlation')
    assert_equal(len(scores), ica.n_components_)

    assert_array_equal(raw_data, raw[:][0])
    assert_array_equal(epochs_data, epochs.get_data())
    assert_array_equal(evoked_data, evoked.data)

    # check score funcs
    for name, func in get_score_funcs().items():

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]
ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5, picks=picks)

ecg_inds, scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
ica.plot_scores(scores, exclude=ecg_inds, title=title % 'ecg', labels='ecg')

show_picks = np.abs(scores).argsort()[::-1][:5]

ica.plot_sources(raw, show_picks, exclude=ecg_inds, title=title % 'ecg')
ica.plot_components(ecg_inds, title=title % 'ecg', colorbar=True)

ecg_inds = ecg_inds[:n_max_ecg]
ica.exclude += ecg_inds

# detect EOG by correlation

eog_inds, scores = ica.find_bads_eog(raw)
ica.plot_scores(scores, exclude=eog_inds, title=title % 'eog', labels='eog')

show_picks = np.abs(scores).argsort()[::-1][:5]

ica.plot_sources(raw, show_picks, exclude=eog_inds, title=title % 'eog')
ica.plot_components(eog_inds, title=title % 'eog', colorbar=True)

eog_inds = eog_inds[:n_max_eog]
ica.exclude += eog_inds

###############################################################################
# 3) Assess component selection and unmixing quality

# estimate average artifact
ecg_evoked = ecg_epochs.average()

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]
event_ids = {"faces": 1, "scrambled": 2}

tmin, tmax = -0.2, 0.6
baseline = None  # no baseline as high-pass is applied
reject = dict(mag=5e-12)

epochs = mne.Epochs(raw, events, event_ids, tmin, tmax,  picks=picks,
                    baseline=baseline, preload=True, reject=reject)

# Fit ICA, find and remove major artifacts
ica = ICA(n_components=0.95, random_state=0).fit(raw, decim=1, reject=reject)

# compute correlation scores, get bad indices sorted by score
eog_epochs = create_eog_epochs(raw, ch_name='MRT31-2908', reject=reject)
eog_inds, eog_scores = ica.find_bads_eog(eog_epochs, ch_name='MRT31-2908')
ica.plot_scores(eog_scores, eog_inds)  # see scores the selection is based on
ica.plot_components(eog_inds)  # view topographic sensitivity of components
ica.exclude += eog_inds[:1]  # we saw the 2nd ECG component looked too dipolar
ica.plot_overlay(eog_epochs.average())  # inspect artifact removal
ica.apply(epochs)  # clean data, default in place

evoked = [epochs[k].average() for k in event_ids]

contrast = combine_evoked(evoked, weights=[-1, 1])  # Faces - scrambled

evoked.append(contrast)

for e in evoked:
    e.plot(ylim=dict(mag=[-400, 400]))

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]
ica = ICA(n_components=0.90, max_pca_components=None)

###############################################################################
# 1) Fit ICA model and identify bad sources

picks = mne.pick_types(raw.info, meg=True, eeg=False, eog=False,
                       stim=False, exclude='bads')

ica.fit(raw, picks=picks, decim=3, reject=dict(mag=4e-12, grad=4000e-13))

# create EOG epochs to improve detection by correlation
picks = mne.pick_types(raw.info, meg=True, eog=True)
eog_epochs = create_eog_epochs(raw, picks=picks)

eog_inds, scores = ica.find_bads_eog(eog_epochs)  # inds sorted!

ica.plot_scores(scores, exclude=eog_inds)  # inspect metrics used

show_picks = np.abs(scores).argsort()[::-1][:5]  # indices of top five scores

# detected artifacts drawn in red (via exclude)
ica.plot_sources(raw, show_picks, exclude=eog_inds, start=0., stop=3.0)
ica.plot_components(eog_inds, colorbar=False)  # show component sensitivites

ica.exclude += eog_inds[:1]  # mark first for exclusion

###############################################################################
# 3) check detection and visualize artifact rejection

# estimate average artifact

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]
def ica_method(raw, picks, plot='n', save ='n'):
    
    
    ###############################################################################
    # 1) Fit ICA model using the FastICA algorithm
    # Other available choices are `infomax` or `extended-infomax`
    # We pass a float value between 0 and 1 to select n_components based on the
    # percentage of variance explained by the PCA components.
    
    ica = ICA(n_components=0.95, method='fastica')
    
    picks = mne.pick_types(raw.info, meg=False, eeg=True, eog=False,
                           stim=False, exclude='bads')
    
    ica.fit(raw, picks=picks, decim=3)
    
    # maximum number of components to reject
    n_max_eog =  1  # here we don't expect horizontal EOG components
    
    ###############################################################################
    # 2) identify bad components by analyzing latent sources.
    
    
    # detect EOG by correlation
    
    eog_inds, scores = ica.find_bads_eog(raw, threshold=2.5)    
    show_picks = np.abs(scores).argsort()[::-1][:5]
    eog_inds = eog_inds[:n_max_eog]
    ica.exclude += eog_inds
    
    ###############################################################################
    # 3) Assess component selection and unmixing quality
    eog_evoked = create_eog_epochs(raw, tmin=-.5, tmax=.5, picks=picks).average()
    
    if plot=='y':
        
        title = 'Sources related to %s artifacts (red)'
        ica.plot_scores(scores, exclude=eog_inds, title=title % 'eog', labels='eog')
        if save=='y':
            pylab.savefig('2.png')

        ica.plot_sources(raw, show_picks, exclude=eog_inds, title=title % 'eog')
        if save=='y':
            pylab.savefig('3.png')

        ica.plot_components(eog_inds, title=title % 'eog', colorbar=True)
        if save=='y':

            pylab.savefig('4.png')

        ica.plot_overlay(raw)  # EOG artifacts remain
        if save=='y':

            pylab.savefig('5.png')

        ica.plot_sources(eog_evoked, exclude=eog_inds)  # plot EOG sources + selection
        if save=='y':

            pylab.savefig('6.png')

        ica.plot_overlay(eog_evoked, exclude=eog_inds)  # plot EOG cleaning
        if save=='y':
            pylab.savefig('7.png')


    ica.apply(raw, exclude=eog_inds)    
    eeg_only_after=raw.pick_types(meg=False, eeg=True)    
    

    return eeg_only_after
    

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]
for trial_type in ['VS']:

    epochs = mne.Epochs(raw, eve_dict[trial_type], id_dict,
                        tmin, tmax, picks=picks, verbose=False,
                        baseline=baseline, reject=reject, preload=True,
                        reject_tmin=rej_tmin, reject_tmax=rej_tmax) # Check rejection settings

    picks = pick_types(epochs.info, meg=True, eog=False)

    # cannot compute EOG covariance?!?!!!
    #baseline_cov = compute_covariance(epochs, tmin=None, tmax=0)
    #ica = ICA(n_components=n_components, max_pca_components=None,max_iter=256, noise_cov=baseline_cov)
    ica = ICA(n_components=n_components, max_pca_components=None,max_iter=256)
    ica.fit(epochs, picks=picks, decim = 5, reject=ica_reject)

    eog_inds, scores = ica.find_bads_eog(epochs, ch_name='EOG001,EOG003')
    allscores = np.vstack((scores[0], scores[1]))
    mscores = np.max(np.abs(allscores), axis=0)
    # now scores is 
    show_picks = mscores.argsort()[::-1][:5]
    eog_inds = list(show_picks[:n_max_eog])
    ica.exclude += eog_inds

    fig = ica.plot_scores(scores, exclude=eog_inds, title=title % 'eog')
    fig.savefig(img_folder + '/ica_eog_scores.png')

    fig = ica.plot_sources(epochs, show_picks, exclude=eog_inds, title=title % 'eog')
    fig.savefig(img_folder + '/ica_eog_sources.png')
    fig = ica.plot_components(show_picks, title=title % 'eog', colorbar=True)
    fig.set_size_inches(12.,8.)
    fig.savefig(img_folder + '/ica_eog_components.png')

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]
# uncomment the code below to test the inteactive mode of plot_components:
# ica.plot_components(picks=range(10), inst=raw)

###############################################################################
# Advanced artifact detection
# ---------------------------
#
# Let's use a more efficient way to find artefacts

eog_average = create_eog_epochs(raw, reject=dict(mag=5e-12, grad=4000e-13),
                                picks=picks_meg).average()

# We simplify things by setting the maximum number of components to reject
n_max_eog = 1  # here we bet on finding the vertical EOG components
eog_epochs = create_eog_epochs(raw, reject=reject)  # get single EOG trials
eog_inds, scores = ica.find_bads_eog(eog_epochs)  # find via correlation

ica.plot_scores(scores, exclude=eog_inds)  # look at r scores of components
# we can see that only one component is highly correlated and that this
# component got detected by our correlation analysis (red).

ica.plot_sources(eog_average, exclude=eog_inds)  # look at source time course

###############################################################################
# We can take a look at the properties of that component, now using the
# data epoched with respect to EOG events.
# We will also use a little bit of smoothing along the trials axis in the
# epochs image:
ica.plot_properties(eog_epochs, picks=eog_inds, psd_args={'fmax': 35.},
                    image_args={'sigma': 1.})

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]
#ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5, picks=picks)

#ecg_inds, scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
#ica.plot_scores(scores, exclude=ecg_inds, title=title % 'ecg', labels='ecg')

#show_picks = np.abs(scores).argsort()[::-1][:5]

#ica.plot_sources(raw, show_picks, exclude=ecg_inds, title=title % 'ecg')
#ica.plot_components(ecg_inds, title=title % 'ecg', colorbar=True)

#ecg_inds = ecg_inds[:n_max_ecg]
#ica.exclude += ecg_inds

# detect EOG by correlation

eog_inds, scores = ica.find_bads_eog(raw, threshold=2.0)
ica.plot_scores(scores, exclude=eog_inds, title=title % 'eog', labels='eog')

show_picks = np.abs(scores).argsort()[::-1][:5]

ica.plot_sources(raw, show_picks, exclude=eog_inds, title=title % 'eog')
ica.plot_components(eog_inds, title=title % 'eog', colorbar=True)

eog_inds = eog_inds[:n_max_eog]
ica.exclude += eog_inds

###############################################################################
# 3) Assess component selection and unmixing quality

# estimate average artifact
#ecg_evoked = ecg_epochs.average()

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]
def runICA(raw,saveRoot,name):

    saveRoot = saveRoot    
    icaList = [] 
    ica = []
    n_max_ecg = 3   # max number of ecg components 
#    n_max_eog_1 = 2 # max number of vert eog comps
#    n_max_eog_2 = 2 # max number of horiz eog comps          
    ecg_source_idx, ecg_scores, ecg_exclude = [], [], []
    eog_source_idx, eog_scores, eog_exclude = [], [], []
    #horiz = 1       # will later be modified to horiz = 0 if no horizontal EOG components are identified                   
    ica = ICA(n_components=0.90,n_pca_components=64,max_pca_components=100,noise_cov=None)
        
    ica.fit(raw)
    #*************
    eog_picks = mne.pick_types(raw.info, meg=False, eeg=False, stim=False, eog=True, ecg=False, emg=False)[0]
    ecg_picks = mne.pick_types(raw.info, meg=False, eeg=False, stim=False, ecg=True, eog=False, emg=False)[0]
    ica_picks = mne.pick_types(raw.info, meg=True, eeg=False, eog=False, ecg=False,
                   stim=False, exclude='bads')
    ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5, picks=ica_picks)
    ecg_evoked = ecg_epochs.average()
    eog_evoked = create_eog_epochs(raw, tmin=-.5, tmax=.5, picks=ica_picks).average()

    ecg_source_idx, ecg_scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
    eog_source_idx, eog_scores = ica.find_bads_eog(raw,ch_name=raw.ch_names[eog_picks].encode('ascii', 'ignore'))
       
    # defining a title-frame for later use
    title = 'Sources related to %s artifacts (red)'

    # extracting number of ica-components and plotting their topographies
    source_idx = range(0, ica.n_components_)
    ica_plot = ica.plot_components(source_idx, ch_type="mag")                                           

    # select ICA sources and reconstruct MEG signals, compute clean ERFs
    # Add detected artefact sources to exclusion list
    # We now add the eog artefacts to the ica.exclusion list
    if not ecg_source_idx:
        print("No ECG components above threshold were identified for subject " + name +
        " - selecting the component with the highest score under threshold")
        ecg_exclude = [np.absolute(ecg_scores).argmax()]
        ecg_source_idx=[np.absolute(ecg_scores).argmax()]
    elif ecg_source_idx:
        ecg_exclude += ecg_source_idx[:n_max_ecg]
    ica.exclude += ecg_exclude

    if not eog_source_idx:
        if np.absolute(eog_scores).any>0.3:
            eog_exclude=[np.absolute(eog_scores).argmax()]
            eog_source_idx=[np.absolute(eog_scores).argmax()]
            print("No EOG components above threshold were identified " + name +
            " - selecting the component with the highest score under threshold above 0.3")
        elif not np.absolute(eog_scores).any>0.3:
            eog_exclude=[]
            print("No EOG components above threshold were identified" + name)
    elif eog_source_idx:
         eog_exclude += eog_source_idx

    ica.exclude += eog_exclude

    print('########## saving')
    if len(eog_exclude) == 0:
        if len(ecg_exclude) == 0:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg_none' + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 1:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg' + map(str, ecg_exclude)[0] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 2:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 3:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '_' + map(str, ecg_exclude)[2] + '.pdf', format = 'pdf')
    elif len(eog_exclude) == 1:
        if len(ecg_exclude) == 0:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg_none' + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 1:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg' + map(str, ecg_exclude)[0] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 2:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 3:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '_' + map(str, ecg_exclude)[2] + '.pdf', format = 'pdf')
    elif len(eog_exclude) == 2:
        if len(ecg_exclude) == 0:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] + '_' + map(str, eog_exclude)[1] +
            '-ecg_none' + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 1:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] + '_' + map(str, eog_exclude)[1] +
            '-ecg' + map(str, ecg_exclude)[0] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 2:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] + '_' + map(str, eog_exclude)[1] +
            '-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 3:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] + '_' + map(str, eog_exclude)[1] +
            '-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '_' + map(str, ecg_exclude)[2] + '.pdf', format = 'pdf')
    
    # plot the scores for the different components highlighting in red that/those related to ECG
    scores_plots_ecg=ica.plot_scores(ecg_scores, exclude=ecg_source_idx, title=title % 'ecg')
    scores_plots_ecg.savefig(saveRoot + name + '_ecg_scores.pdf', format = 'pdf')
    scores_plots_eog=ica.plot_scores(eog_scores, exclude=eog_source_idx, title=title % 'eog')
#.........这里部分代码省略.........

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]
def compute_ica(subject, data_folder):
    """Function will compute ICA on raw and apply the ICA.

    Parameters
    ----------
    subject : string
        the subject id to be loaded
    """
    raw = mne.io.Raw(data_folder + "%s_bp-raw.fif" % subject, preload=True)
    raw.set_montage = montage
    raw.apply_proj()
    # raw.resample(512, n_jobs=2)

    # ICA Part
    ica = ICA(n_components=None,
              max_pca_components=40,
              method='fastica',
              max_iter=256)

    picks = mne.pick_types(
        raw.info, meg=False, eeg=True, stim=False, exclude='bads')

    ica.fit(raw, picks=picks, decim=decim, reject=reject)

    # maximum number of components to reject
    n_max_eog = 1

    ##########################################################################
    # 2) identify bad components by analyzing latent sources.
    title = 'Sources related to %s artifacts (red) for sub: %s'
    #
    # # generate ECG epochs use detection via phase statistics
    # ecg_epochs = create_ecg_epochs(raw, ch_name="Ext4",
    #                                tmin=-.5, tmax=.5, picks=picks)
    # n_ecg_epochs_found = len(ecg_epochs.events)
    # sel_ecg_epochs = np.arange(0, n_ecg_epochs_found, 10)
    # ecg_epochs = ecg_epochs[sel_ecg_epochs]
    #
    # ecg_inds, scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
    # fig = ica.plot_scores(scores, exclude=ecg_inds,
    #                       title=title % ('ecg', subject))
    # fig.savefig(data_folder + "pics/%s_ecg_scores.png" % subject)
    #
    # if ecg_inds:
    #     show_picks = np.abs(scores).argsort()[::-1][:5]
    #
    #     fig = ica.plot_sources(raw, show_picks, exclude=ecg_inds,
    #                            title=title % ('ecg', subject), show=False)
    #     fig.savefig(data_folder + "pics/%s_ecg_sources.png" % subject)
    #     fig = ica.plot_components(ecg_inds, title=title % ('ecg', subject),
    #                               colorbar=True)
    #     fig.savefig(data_folder + "pics/%s_ecg_component.png" % subject)
    #
    #     ecg_inds = ecg_inds[:n_max_ecg]
    #     ica.exclude += ecg_inds
    #
    # # estimate average artifact
    # ecg_evoked = ecg_epochs.average()
    # del ecg_epochs
    #
    # # plot ECG sources + selection
    # fig = ica.plot_sources(ecg_evoked, exclude=ecg_inds)
    # fig.savefig(data_folder + "pics/%s_ecg_sources_ave.png" % subject)
    #
    # # plot ECG cleaning
    # ica.plot_overlay(ecg_evoked, exclude=ecg_inds)
    # fig.savefig(data_folder + "pics/%s_ecg_sources_clean_ave.png" % subject)

    # DETECT EOG BY CORRELATION
    # HORIZONTAL EOG
    eog_epochs = create_eog_epochs(raw, ch_name="EXG4")
    eog_indices, scores = ica.find_bads_eog(raw, ch_name="EXG4")
    fig = ica.plot_scores(
        scores, exclude=eog_indices, title=title % ('eog', subject))
    fig.savefig(data_folder + "pics/%s_eog_scores.png" % subject)

    fig = ica.plot_components(
        eog_indices, title=title % ('eog', subject), colorbar=True)
    fig.savefig(data_folder + "pics/%s_eog_component.png" % subject)

    eog_indices = eog_indices[:n_max_eog]
    ica.exclude += eog_indices

    del eog_epochs

    ##########################################################################
    # Apply the solution to Raw, Epochs or Evoked like this:
    raw_ica = ica.apply(raw)
    ica.save(data_folder + "%s-ica.fif" % subject)  # save ICA componenets
    # Save raw with ICA removed
    raw_ica.save(data_folder + "%s_bp_ica-raw.fif" % subject, overwrite=True)
    plt.close("all")

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]
def runICA(raw,saveRoot,name):

    saveRoot = saveRoot    
    icaList = [] 
    ica = []
    n_max_ecg = 3   # max number of ecg components 
#    n_max_eog_1 = 2 # max number of vert eog comps
#    n_max_eog_2 = 2 # max number of horiz eog comps          
    ecg_source_idx, ecg_scores, ecg_exclude = [], [], []
    eog_source_idx, eog_scores, eog_exclude = [], [], []
    #horiz = 1       # will later be modified to horiz = 0 if no horizontal EOG components are identified                   
    ica = ICA(n_components=0.90,n_pca_components=64,max_pca_components=100,noise_cov=None)
    
    fit_picks = mne.pick_types(raw.info, meg=True, eeg=True, eog=False, ecg=False,
                   stim=False, exclude='bads')    
    ica.fit(raw, picks=fit_picks)
    #ica.fit(raw)
    #*************
    eog_picks = mne.pick_types(raw.info, meg=False, eeg=False, stim=False, eog=True, ecg=False, emg=False)[0]
    ecg_picks = mne.pick_types(raw.info, meg=False, eeg=False, stim=False, ecg=True, eog=False, emg=False)[0]
    ica_picks = mne.pick_types(raw.info, meg=True, eeg=True, eog=False, ecg=False,
                   stim=False, exclude='bads')
    ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5, picks=ica_picks)
    ecg_evoked = ecg_epochs.average()
    eog_evoked = create_eog_epochs(raw, tmin=-.5, tmax=.5, picks=ica_picks).average()

    ecg_source_idx, ecg_scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
    eog_source_idx, eog_scores = ica.find_bads_eog(raw,ch_name=raw.ch_names[eog_picks].encode('ascii', 'ignore'))
    #eog_source_idx_2, eog_scores_2 = ica.find_bads_eog(raw,ch_name='EOG002')
    #if not eog_source_idx_2:
    #    horiz = 0
    
    #show_picks = np.abs(scores).argsort()[::-1][:5]
    #ica.plot_sources(raw, show_picks, exclude=ecg_inds, title=title % 'ecg')
    
        
    # defining a title-frame for later use
    title = 'Sources related to %s artifacts (red)'
  

    # extracting number of ica-components and plotting their topographies
    source_idx = range(0, ica.n_components_)
    #ica_plot = ica.plot_components(source_idx, ch_type="mag")
    ica_plot = ica.plot_components(source_idx)
                                          
    #ica_plot = ica.plot_components(source_idx)

    # select ICA sources and reconstruct MEG signals, compute clean ERFs
    # Add detected artefact sources to exclusion list
    # We now add the eog artefacts to the ica.exclusion list
    if not ecg_source_idx:
        print("No ECG components above threshold were identified for subject " + name +
        " - selecting the component with the highest score under threshold")
        ecg_exclude = [np.absolute(ecg_scores).argmax()]
        ecg_source_idx=[np.absolute(ecg_scores).argmax()]
    elif ecg_source_idx:
        ecg_exclude += ecg_source_idx[:n_max_ecg]
    ica.exclude += ecg_exclude

    if not eog_source_idx:
        if np.absolute(eog_scores).any>0.3:
            eog_exclude=[np.absolute(eog_scores).argmax()]
            eog_source_idx=[np.absolute(eog_scores).argmax()]
            print("No EOG components above threshold were identified " + name +
            " - selecting the component with the highest score under threshold above 0.3")
        elif not np.absolute(eog_scores).any>0.3:
            eog_exclude=[]
            print("No EOG components above threshold were identified" + name)
    elif eog_source_idx:
         eog_exclude += eog_source_idx

    ica.exclude += eog_exclude

    print('########## saving')
    if len(eog_exclude) == 0:
        if len(ecg_exclude) == 0:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg_none' + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 1:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg' + map(str, ecg_exclude)[0] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 2:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 3:
            ica_plot.savefig(saveRoot + name + '_comps_eog_none-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '_' + map(str, ecg_exclude)[2] + '.pdf', format = 'pdf')
    elif len(eog_exclude) == 1:
        if len(ecg_exclude) == 0:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg_none' + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 1:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg' + map(str, ecg_exclude)[0] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 2:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '.pdf', format = 'pdf')
        elif len(ecg_exclude) == 3:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] +
            '-ecg' + map(str, ecg_exclude)[0] + '_' + map(str, ecg_exclude)[1] + '_' + map(str, ecg_exclude)[2] + '.pdf', format = 'pdf')
    elif len(eog_exclude) == 2:
        if len(ecg_exclude) == 0:
            ica_plot.savefig(saveRoot + name + '_comps_eog' + map(str, eog_exclude)[0] + '_' + map(str, eog_exclude)[1] +
            '-ecg_none' + '.pdf', format = 'pdf')
#.........这里部分代码省略.........

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]
def compute_ica(subject):
    """Function will compute ICA on raw and apply the ICA.

    params:
    subject : str
        the subject id to be loaded
    """
    raw = Raw(save_folder + "%s_filtered_data_mc_raw_tsss.fif" % subject,
              preload=True)

    # ICA Part
    ica = ICA(n_components=0.95, method='fastica', max_iter=256)

    picks = mne.pick_types(raw.info, meg=True, eeg=True,
                           stim=False, exclude='bads')

    ica.fit(raw, picks=picks, decim=decim, reject=reject)

    # maximum number of components to reject
    n_max_ecg, n_max_eog = 3, 1

    ##########################################################################
    # 2) identify bad components by analyzing latent sources.
    title = 'Sources related to %s artifacts (red) for sub: %s'

    # generate ECG epochs use detection via phase statistics
    ecg_epochs = create_ecg_epochs(raw, ch_name="ECG002",
                                   tmin=-.5, tmax=.5, picks=picks)
    n_ecg_epochs_found = len(ecg_epochs.events)
    sel_ecg_epochs = np.arange(0, n_ecg_epochs_found, 10)
    ecg_epochs = ecg_epochs[sel_ecg_epochs]

    ecg_inds, scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
    fig = ica.plot_scores(scores, exclude=ecg_inds,
                          title=title % ('ecg', subject))
    fig.savefig(save_folder + "pics/%s_ecg_scores.png" % subject)

    if ecg_inds:
        show_picks = np.abs(scores).argsort()[::-1][:5]

        fig = ica.plot_sources(raw, show_picks, exclude=ecg_inds,
                               title=title % ('ecg', subject), show=False)
        fig.savefig(save_folder + "pics/%s_ecg_sources.png" % subject)
        fig = ica.plot_components(ecg_inds, title=title % ('ecg', subject),
                                  colorbar=True)
        fig.savefig(save_folder + "pics/%s_ecg_component.png" % subject)

        ecg_inds = ecg_inds[:n_max_ecg]
        ica.exclude += ecg_inds

    # estimate average artifact
    ecg_evoked = ecg_epochs.average()
    del ecg_epochs

    # plot ECG sources + selection
    fig = ica.plot_sources(ecg_evoked, exclude=ecg_inds)
    fig.savefig(save_folder + "pics/%s_ecg_sources_ave.png" % subject)

    # plot ECG cleaning
    ica.plot_overlay(ecg_evoked, exclude=ecg_inds)
    fig.savefig(save_folder + "pics/%s_ecg_sources_clean_ave.png" % subject)

    # DETECT EOG BY CORRELATION
    # HORIZONTAL EOG
    eog_epochs = create_eog_epochs(raw, ch_name="EOG001")
    eog_inds, scores = ica.find_bads_eog(raw)
    fig = ica.plot_scores(scores, exclude=eog_inds,
                          title=title % ('eog', subject))
    fig.savefig(save_folder + "pics/%s_eog_scores.png" % subject)

    fig = ica.plot_components(eog_inds, title=title % ('eog', subject),
                              colorbar=True)
    fig.savefig(save_folder + "pics/%s_eog_component.png" % subject)

    eog_inds = eog_inds[:n_max_eog]
    ica.exclude += eog_inds

    del eog_epochs

    ##########################################################################
    # Apply the solution to Raw, Epochs or Evoked like this:
    raw_ica = ica.apply(raw, copy=False)
    ica.save(save_folder + "%s-ica.fif" % subject)  # save ICA componenets
    # Save raw with ICA removed
    raw_ica.save(save_folder + "%s_filtered_ica_mc_raw_tsss.fif" % subject,
                 overwrite=True)
    plt.close("all")

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]

#.........这里部分代码省略.........
    # check scrore funcs
    for name, func in score_funcs.items():
        if name in score_funcs_unsuited:
            continue
        scores = ica.score_sources(raw, target="EOG 061", score_func=func, start=0, stop=10)
        assert_true(ica.n_components_ == len(scores))

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

    params = []
    params += [(None, -1, slice(2), [0, 1])]  # varicance, kurtosis idx params
    params += [(None, "MEG 1531")]  # ECG / EOG channel params
    for idx, ch_name in product(*params):
        ica.detect_artifacts(
            raw,
            start_find=0,
            stop_find=50,
            ecg_ch=ch_name,
            eog_ch=ch_name,
            skew_criterion=idx,
            var_criterion=idx,
            kurt_criterion=idx,
        )

    idx, scores = ica.find_bads_ecg(raw, method="ctps")
    assert_equal(len(scores), ica.n_components_)
    idx, scores = ica.find_bads_ecg(raw, method="correlation")
    assert_equal(len(scores), ica.n_components_)
    idx, scores = ica.find_bads_ecg(epochs, method="ctps")
    assert_equal(len(scores), ica.n_components_)
    assert_raises(ValueError, ica.find_bads_ecg, epochs.average(), method="ctps")
    assert_raises(ValueError, ica.find_bads_ecg, raw, method="crazy-coupling")

    idx, scores = ica.find_bads_eog(raw)
    assert_equal(len(scores), ica.n_components_)
    raw.info["chs"][raw.ch_names.index("EOG 061") - 1]["kind"] = 202
    idx, scores = ica.find_bads_eog(raw)
    assert_true(isinstance(scores, list))
    assert_equal(len(scores[0]), ica.n_components_)

    # check score funcs
    for name, func in score_funcs.items():
        if name in score_funcs_unsuited:
            continue
        scores = ica.score_sources(epochs_eog, target="EOG 061", score_func=func)
        assert_true(ica.n_components_ == len(scores))

    # check univariate stats
    scores = ica.score_sources(epochs, score_func=stats.skew)

    # check exception handling
    assert_raises(ValueError, ica.score_sources, epochs, target=np.arange(1))

    # ecg functionality
    ecg_scores = ica.score_sources(raw, target="MEG 1531", score_func="pearsonr")

    with warnings.catch_warnings(record=True):  # filter attenuation warning
        ecg_events = ica_find_ecg_events(raw, sources[np.abs(ecg_scores).argmax()])

    assert_true(ecg_events.ndim == 2)

    # eog functionality
    eog_scores = ica.score_sources(raw, target="EOG 061", score_func="pearsonr")
    with warnings.catch_warnings(record=True):  # filter attenuation warning
        eog_events = ica_find_eog_events(raw, sources[np.abs(eog_scores).argmax()])

    assert_true(eog_events.ndim == 2)

    # Test ica fiff export
    ica_raw = ica.get_sources(raw, start=0, stop=100)
    assert_true(ica_raw.last_samp - ica_raw.first_samp == 100)
    assert_true(len(ica_raw._filenames) == 0)  # API consistency
    ica_chans = [ch for ch in ica_raw.ch_names if "ICA" in ch]
    assert_true(ica.n_components_ == len(ica_chans))
    test_ica_fname = op.join(op.abspath(op.curdir), "test-ica_raw.fif")
    ica.n_components = np.int32(ica.n_components)
    ica_raw.save(test_ica_fname, overwrite=True)
    ica_raw2 = io.Raw(test_ica_fname, preload=True)
    assert_allclose(ica_raw._data, ica_raw2._data, rtol=1e-5, atol=1e-4)
    ica_raw2.close()
    os.remove(test_ica_fname)

    # Test ica epochs export
    ica_epochs = ica.get_sources(epochs)
    assert_true(ica_epochs.events.shape == epochs.events.shape)
    ica_chans = [ch for ch in ica_epochs.ch_names if "ICA" in ch]
    assert_true(ica.n_components_ == len(ica_chans))
    assert_true(ica.n_components_ == ica_epochs.get_data().shape[1])
    assert_true(ica_epochs.raw is None)
    assert_true(ica_epochs.preload is True)

    # test float n pca components
    ica.pca_explained_variance_ = np.array([0.2] * 5)
    ica.n_components_ = 0
    for ncomps, expected in [[0.3, 1], [0.9, 4], [1, 1]]:
        ncomps_ = _check_n_pca_components(ica, ncomps)
        assert_true(ncomps_ == expected)

# 需要导入模块: from mne.preprocessing import ICA [as 别名]
# 或者: from mne.preprocessing.ICA import find_bads_eog [as 别名]
def preprocess_ICA_fif_to_ts(fif_file, ECG_ch_name, EoG_ch_name, l_freq, h_freq):
    # ------------------------ Import stuff ------------------------ #
    import os
    import mne
    import sys
    from mne.io import Raw
    from mne.preprocessing import ICA
    from mne.preprocessing import create_ecg_epochs, create_eog_epochs
    from nipype.utils.filemanip import split_filename as split_f
    from reportGen import generateReport
    import pickle

    subj_path, basename, ext = split_f(fif_file)
    # -------------------- Delete later ------------------- #
    subj_name = subj_path[-5:]
    results_dir = subj_path[:-6]
    # results_dir += '2016'
    subj_path = results_dir + '/' + subj_name
    if not os.path.exists(subj_path):
        try:
            os.makedirs(subj_path)
        except OSError:
            sys.stderr.write('ica_preproc: problem creating directory: ' + subj_path)
    ########################################################
    # Read raw
    #   If None the compensation in the data is not modified. If set to n, e.g. 3, apply
    #   gradient compensation of grade n as for CTF systems (compensation=3)
    print(fif_file)
    print(EoG_ch_name)
    #  ----------------------------- end Import stuff ----------------- #
    # EoG_ch_name = "EOG061, EOG062"

    # ------------- Load raw ------------- #
    raw = Raw(fif_file, preload=True)
    # select sensors
    select_sensors = mne.pick_types(raw.info, meg=True, ref_meg=False, exclude='bads')
    picks_meeg = mne.pick_types(raw.info, meg=True, eeg=True, exclude='bads')

    # filtering
    raw.filter(l_freq=l_freq, h_freq=h_freq, picks=picks_meeg, method='iir', n_jobs=1)

    # if ECG_ch_name == 'EMG063':
    if ECG_ch_name in raw.info['ch_names']:
        raw.set_channel_types({ECG_ch_name: 'ecg'})  # Without this files with ECG_ch_name = 'EMG063' fail
        # ECG_ch_name = 'ECG063'
    if EoG_ch_name == 'EMG065,EMG066,EMG067,EMG068':   # Because ica.find_bads_eog... can process max 2 EoG channels
        EoG_ch_name = 'EMG065,EMG067'                 # it won't fail if I specify 4 channels, but it'll use only first
                                                      # EMG065 and EMG066 are for vertical eye movements and
                                                      # EMG067 and EMG068 are for horizontal

    # print rnk
    rnk = 'N/A'
    # 1) Fit ICA model using the FastICA algorithm
    # Other available choices are `infomax` or `extended-infomax`
    # We pass a float value between 0 and 1 to select n_components based on the
    # percentage of variance explained by the PCA components.
    reject = dict(mag=10e-12, grad=10000e-13)
    flat = dict(mag=0.1e-12, grad=1e-13)
    # check if we have an ICA, if yes, we load it
    ica_filename = os.path.join(subj_path, basename + "-ica.fif")
    raw_ica_filename = os.path.join(subj_path, basename + "_ica_raw.fif")
    info_filename = os.path.join(subj_path, basename + "_info.pickle")
    # if os.path.exists(ica_filename) == False:
    ica = ICA(n_components=0.99, method='fastica')  # , max_iter=500
    ica.fit(raw, picks=select_sensors, reject=reject, flat=flat)  # decim = 3,
    # has_ICA = False
    # else:
    #     has_ICA = True
    #     ica = read_ica(ica_filename)
    #     ica.exclude = []
    # ica.labels_ = dict() # to avoid bug; Otherwise it'll throw an exception

    ica_sources_filename = subj_path + '/' + basename + '_ica_timecourse.fif'

    # if not os.path.isfile(ica_sources_filename):
    icaSrc = ica.get_sources(raw, add_channels=None, start=None, stop=None)
    icaSrc.save(ica_sources_filename, picks=None, tmin=0, tmax=None, buffer_size_sec=10,
                drop_small_buffer=False, proj=False, fmt='single', overwrite=True, split_size='2GB', verbose=None)
    icaSrc = None
    # if has_ICA == False:
    # ica.save(ica_filename)
    # return
    # 2) identify bad components by analyzing latent sources.
    # generate ECG epochs use detection via phase statistics

    # check if ECG_ch_name is in the raw channels
    # import ipdb; ipdb.set_trace()
    if ECG_ch_name in raw.info['ch_names']:
        ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5, picks=select_sensors, ch_name=ECG_ch_name)
    # if not  a synthetic ECG channel is created from cross channel average
    else:
        ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5, picks=select_sensors)
    # ICA for ECG artifact
    # threshold=0.25 come defualt
    ecg_inds, ecg_scores = ica.find_bads_ecg(ecg_epochs, method='ctps', threshold=0.25)
    # if len(ecg_inds) > 0:
    ecg_evoked = ecg_epochs.average()
    ecg_epochs = None    # ecg_epochs use too much memory
    n_max_ecg = 3
    ecg_inds = ecg_inds[:n_max_ecg]
#.........这里部分代码省略.........