Python preprocessing.create_eog_epochs函数代码示例

def test_plot_ica_overlay():
    """Test plotting of ICA cleaning."""
    import matplotlib.pyplot as plt
    raw = _get_raw(preload=True)
    picks = _get_picks(raw)
    ica = ICA(noise_cov=read_cov(cov_fname), n_components=2,
              max_pca_components=3, n_pca_components=3)
    # can't use info.normalize_proj here because of how and when ICA and Epochs
    # objects do picking of Raw data
    with pytest.warns(RuntimeWarning, match='projection'):
        ica.fit(raw, picks=picks)
    # don't test raw, needs preload ...
    with pytest.warns(RuntimeWarning, match='projection'):
        ecg_epochs = create_ecg_epochs(raw, picks=picks)
    ica.plot_overlay(ecg_epochs.average())
    with pytest.warns(RuntimeWarning, match='projection'):
        eog_epochs = create_eog_epochs(raw, picks=picks)
    ica.plot_overlay(eog_epochs.average())
    pytest.raises(TypeError, ica.plot_overlay, raw[:2, :3][0])
    ica.plot_overlay(raw)
    plt.close('all')

    # smoke test for CTF
    raw = read_raw_fif(raw_ctf_fname)
    raw.apply_gradient_compensation(3)
    picks = pick_types(raw.info, meg=True, ref_meg=False)
    ica = ICA(n_components=2, max_pca_components=3, n_pca_components=3)
    ica.fit(raw, picks=picks)
    with pytest.warns(RuntimeWarning, match='longer than'):
        ecg_epochs = create_ecg_epochs(raw)
    ica.plot_overlay(ecg_epochs.average())
    plt.close('all')

def test_plot_ica_overlay():
    """Test plotting of ICA cleaning
    """
    raw = _get_raw()
    picks = _get_picks(raw)
    ica_picks = pick_types(raw.info, meg=True, eeg=False, stim=False, ecg=False, eog=False, exclude="bads")
    ica = ICA(noise_cov=read_cov(cov_fname), n_components=2, max_pca_components=3, n_pca_components=3)
    ica.fit(raw, picks=ica_picks)
    # don't test raw, needs preload ...
    ecg_epochs = create_ecg_epochs(raw, picks=picks)
    ica.plot_overlay(ecg_epochs.average())
    eog_epochs = create_eog_epochs(raw, picks=picks)
    ica.plot_overlay(eog_epochs.average())
    assert_raises(ValueError, ica.plot_overlay, raw[:2, :3][0])
    plt.close("all")

def test_plot_ica_overlay():
    """Test plotting of ICA cleaning
    """
    raw = _get_raw()
    picks = _get_picks(raw)
    ica = ICA(noise_cov=read_cov(cov_fname), n_components=2,
              max_pca_components=3, n_pca_components=3)
    ica.fit(raw, picks=picks)
    # don't test raw, needs preload ...
    ecg_epochs = create_ecg_epochs(raw, picks=picks)
    ica.plot_overlay(ecg_epochs.average())
    eog_epochs = create_eog_epochs(raw, picks=picks)
    ica.plot_overlay(eog_epochs.average())
    assert_raises(ValueError, ica.plot_overlay, raw[:2, :3][0])
    plt.close('all')

def test_plot_ica_overlay():
    """Test plotting of ICA cleaning."""
    import matplotlib.pyplot as plt
    raw = _get_raw(preload=True)
    picks = _get_picks(raw)
    ica = ICA(noise_cov=read_cov(cov_fname), n_components=2,
              max_pca_components=3, n_pca_components=3)
    # can't use info.normalize_proj here because of how and when ICA and Epochs
    # objects do picking of Raw data
    with warnings.catch_warnings(record=True):  # bad proj
        ica.fit(raw, picks=picks)
    # don't test raw, needs preload ...
    with warnings.catch_warnings(record=True):  # bad proj
        ecg_epochs = create_ecg_epochs(raw, picks=picks)
    ica.plot_overlay(ecg_epochs.average())
    with warnings.catch_warnings(record=True):  # bad proj
        eog_epochs = create_eog_epochs(raw, picks=picks)
    ica.plot_overlay(eog_epochs.average())
    assert_raises(ValueError, ica.plot_overlay, raw[:2, :3][0])
    ica.plot_overlay(raw)
    plt.close('all')

mne.viz.plot_events(events, raw.info['sfreq'])

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]))

# 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.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

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
    

        exclude='bads')

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

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

    ##########################################################################
    # 2) identify bad components by analyzing latent sources.

    # DETECT EOG BY CORRELATION
    # HORIZONTAL EOG
    title = "ICA: %s for %s"

    eog_epochs = create_eog_epochs(raw, ch_name="EOG002")  # TODO: check EOG
    eog_average = eog_epochs.average()
    # channel name
    eog_inds, scores = ica.find_bads_eog(raw)

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

    if eog_inds:
        fig = ica.plot_scores(
            scores, exclude=eog_inds, title=title % ('eog', subject))
        fig.savefig(ica_folder + "plots/%s_%s_eog_scores_2.png" % (subject,
                                                                   condition))
        fig = ica.plot_sources(eog_average, exclude=eog_inds)
        fig.savefig(ica_folder + "plots/%s_%s_eog_source_2.png" % (subject,
                                                                   condition))

#.........这里部分代码省略.........

        fig = ica.plot_scores(scores, exclude=ecg_inds, labels='ecg',
                              title='', show=show)

        report.add_figs_to_section(fig, 'scores ({})'.format(subject),
                                   section=comment + 'ECG',
                                   scale=img_scale)

        current_exclude = [e for e in ica.exclude]  # issue #2608 MNE
        fig = ica.plot_sources(raw, ecg_inds, exclude=ecg_inds,
                               title=title % ('components', 'ecg'), show=show)

        report.add_figs_to_section(fig, 'sources ({})'.format(subject),
                                   section=comment + 'ECG',
                                   scale=img_scale)
        ica.exclude = current_exclude

        fig = ica.plot_components(ecg_inds, ch_type=topo_ch_type,
                                  title='', colorbar=True, show=show)
        report.add_figs_to_section(fig, title % ('sources', 'ecg'),
                                   section=comment + 'ECG', scale=img_scale)
        ica.exclude = current_exclude

        ecg_inds = ecg_inds[:n_max_ecg]
        ica.exclude += ecg_inds
        fig = ica.plot_sources(ecg_evoked, exclude=ecg_inds, show=show)
        report.add_figs_to_section(fig, 'evoked sources ({})'.format(subject),
                                   section=comment + 'ECG',
                                   scale=img_scale)

        fig = ica.plot_overlay(ecg_evoked, exclude=ecg_inds, show=show)
        report.add_figs_to_section(fig,
                                   'rejection overlay ({})'.format(subject),
                                   section=comment + 'ECG',
                                   scale=img_scale)

    # detect EOG by correlation
    picks_eog = np.concatenate(
        [picks_, pick_types(raw.info, meg=False, eeg=False, ecg=False,
                            eog=True)])

    eog_epochs = create_eog_epochs(raw, tmin=eog_tmin, tmax=eog_tmax,
                                   picks=picks_eog, reject=reject_)
    artifact_stats['eog_n_events'] = len(eog_epochs.events)
    artifact_stats['eog_n_used'] = artifact_stats['eog_n_events']
    eog_ave = eog_epochs.average()
    report.add_figs_to_section(eog_ave.plot(), 'EOG-used', 'artifacts')
    _put_artifact_range(artifact_stats, eog_ave, kind='eog')

    eog_inds = None
    if len(eog_epochs.events) > 0:
        eog_inds, scores = ica.find_bads_eog(eog_epochs)

    if eog_inds is not None and len(eog_epochs.events) > 0:
        fig = ica.plot_scores(scores, exclude=eog_inds, labels='eog',
                              show=show, title='')
        report.add_figs_to_section(fig, 'scores ({})'.format(subject),
                                   section=comment + 'EOG',
                                   scale=img_scale)

        current_exclude = [e for e in ica.exclude]  # issue #2608 MNE
        fig = ica.plot_sources(raw, eog_inds, exclude=ecg_inds,
                               title=title % ('sources', 'eog'), show=show)
        report.add_figs_to_section(fig, 'sources', section=comment + 'EOG',
                                   scale=img_scale)
        ica.exclude = current_exclude

        fig = ica.plot_components(eog_inds, ch_type=topo_ch_type,
                                  title='', colorbar=True, show=show)
        report.add_figs_to_section(fig, title % ('components', 'eog'),
                                   section=comment + 'EOG', scale=img_scale)
        ica.exclude = current_exclude

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

        eog_evoked = eog_epochs.average()
        fig = ica.plot_sources(eog_evoked, exclude=eog_inds, show=show)
        report.add_figs_to_section(
            fig, 'evoked sources ({})'.format(subject),
            section=comment + 'EOG', scale=img_scale)

        fig = ica.plot_overlay(eog_evoked, exclude=eog_inds, show=show)
        report.add_figs_to_section(
            fig, 'rejection overlay({})'.format(subject),
            section=comment + 'EOG', scale=img_scale)
    else:
        del eog_epochs

    # check the amplitudes do not change
    if len(ica.exclude) > 0:
        fig = ica.plot_overlay(raw, show=show)  # EOG artifacts remain
        html = _render_components_table(ica)
        report.add_htmls_to_section(
            html, captions='excluded components',
            section='ICA rejection summary (%s)' % ch_type)
        report.add_figs_to_section(
            fig, 'rejection overlay({})'.format(subject),
            section=comment + 'RAW', scale=img_scale)
    return ica, dict(html=report, stats=artifact_stats)

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')
#.........这里部分代码省略.........

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")

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')
#.........这里部分代码省略.........

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")

        icacomps.exclude=[]
        if not icacomps.exclude:
            print('Old components copied. Exclude field cleared')    
    
    raw = mne.io.Raw(rawRoot+name+'.fif', preload=True)
    ecg_picks = mne.pick_types(raw.info, meg=False, eeg=False, eog=False, ecg=True,
                   stim=False, exclude='bads')[0]
    eog_picks = mne.pick_types(raw.info, meg=False, eeg=False, ecg=False, eog=True,
                   stim=False, exclude='bads')[0]
    meg_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=meg_picks, verbose=False)
                                   #ch_name=raw.ch_names[ecg_picks].encode('UTF8'))
    ecg_evoked = ecg_epochs.average()
    eog_evoked = create_eog_epochs(raw, tmin=-.5, tmax=.5,picks=meg_picks,
                           ch_name=raw.ch_names[eog_picks].encode('UTF8'), verbose=False).average()


    # ica topos
    source_idx = range(0, icacomps.n_components_)
    ica_plot = icacomps.plot_components(source_idx, ch_type="mag") 
    plt.waitforbuttonpress(1)
    
    title = 'Sources related to %s artifacts (red)'
    
    #ask for comps ECG
    prompt = '> '
    ecg_done = 'N'
    eog_done = 'N'
    
    while ecg_done.strip() != 'Y' and ecg_done.strip() != 'y':

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]
#.........这里部分代码省略.........