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_find_ecg():
    """Test find ECG peaks."""
    # Test if ECG analysis will work on data that is not preloaded
    raw = read_raw_fif(raw_fname, preload=False)

    # once with mag-trick
    # once with characteristic channel
    for ch_name in ['MEG 1531', None]:
        events, ch_ECG, average_pulse, ecg = find_ecg_events(
            raw, event_id=999, ch_name=ch_name, return_ecg=True)
        assert raw.n_times == ecg.shape[-1]
        n_events = len(events)
        _, times = raw[0, :]
        assert 55 < average_pulse < 60

    picks = pick_types(
        raw.info, meg='grad', eeg=False, stim=False,
        eog=False, ecg=True, emg=False, ref_meg=False,
        exclude='bads')

    # There should be no ECG channels, or else preloading will not be
    # tested
    assert 'ecg' not in raw

    ecg_epochs = create_ecg_epochs(raw, picks=picks, keep_ecg=True)
    assert len(ecg_epochs.events) == n_events
    assert 'ECG-SYN' not in raw.ch_names
    assert 'ECG-SYN' in ecg_epochs.ch_names

    picks = pick_types(
        ecg_epochs.info, meg=False, eeg=False, stim=False,
        eog=False, ecg=True, emg=False, ref_meg=False,
        exclude='bads')
    assert len(picks) == 1

    ecg_epochs = create_ecg_epochs(raw, ch_name='MEG 2641')
    assert 'MEG 2641' in ecg_epochs.ch_names

    # test with user provided ecg channel
    raw.info['projs'] = list()
    with pytest.warns(RuntimeWarning, match='unit for channel'):
        raw.set_channel_types({'MEG 2641': 'ecg'})
    create_ecg_epochs(raw)

    raw.load_data().pick_types()  # remove ECG
    ecg_epochs = create_ecg_epochs(raw, keep_ecg=False)
    assert len(ecg_epochs.events) == n_events
    assert 'ECG-SYN' not in raw.ch_names
    assert 'ECG-SYN' not in ecg_epochs.ch_names

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

#
# License: BSD (3-clause)

import mne
from mne.io import Raw
from mne.preprocessing import ICA, create_ecg_epochs
from mne.datasets import sample

print(__doc__)

###############################################################################
# Fit ICA model using the FastICA algorithm, detect and inspect components

data_path = sample.data_path()
raw_fname = data_path + "/MEG/sample/sample_audvis_filt-0-40_raw.fif"

raw = Raw(raw_fname, preload=True)
raw.filter(1, 30, method="iir")
raw.pick_types(meg=True, eeg=False, exclude="bads", stim=True)

# longer + more epochs for more artifact exposure
events = mne.find_events(raw, stim_channel="STI 014")
epochs = mne.Epochs(raw, events, event_id=None, tmin=-0.2, tmax=0.5)

ica = ICA(n_components=0.95, method="fastica").fit(epochs)

ecg_epochs = create_ecg_epochs(raw, tmin=-0.5, tmax=0.5)
ecg_inds, scores = ica.find_bads_ecg(ecg_epochs)

ica.plot_components(ecg_inds)

def preprocess_ICA_fif_to_ts(fif_file, ECG_ch_name, EoG_ch_name, l_freq, h_freq, down_sfreq, variance, is_sensor_space, data_type):
    import os
    import numpy as np

    import mne
    from mne.io import Raw
    from mne.preprocessing import ICA, read_ica
    from mne.preprocessing import create_ecg_epochs, create_eog_epochs
    from mne.report import Report

    from nipype.utils.filemanip import split_filename as split_f

    report = Report()

    subj_path, basename, ext = split_f(fif_file)
    (data_path, sbj_name) = os.path.split(subj_path)
    print data_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)
    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')

    # save electrode locations
    sens_loc = [raw.info['chs'][i]['loc'][:3] for i in select_sensors]
    sens_loc = np.array(sens_loc)

    channel_coords_file = os.path.abspath("correct_channel_coords.txt")
    print '*** ' + channel_coords_file + '***'
    np.savetxt(channel_coords_file, sens_loc, fmt='%s')

    # save electrode names
    sens_names = np.array([raw.ch_names[pos] for pos in select_sensors],dtype = "str")

    # AP 21032016 
#    channel_names_file = os.path.join(data_path, "correct_channel_names.txt") 
    channel_names_file = os.path.abspath("correct_channel_names.txt")
    np.savetxt(channel_names_file,sens_names , fmt = '%s')
 
    ### filtering + downsampling
    raw.filter(l_freq=l_freq, h_freq=h_freq, picks=picks_meeg,
               method='iir', n_jobs=8)
#    raw.filter(l_freq = l_freq, h_freq = h_freq, picks = picks_meeg,
#               method='iir')
#    raw.resample(sfreq=down_sfreq, npad=0)

    ### 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_title = 'Sources related to %s artifacts (red)'
    is_show = False # visualization
    reject = dict(mag=4e-12, grad=4000e-13)

    # check if we have an ICA, if yes, we load it
    ica_filename = os.path.join(subj_path,basename + "-ica.fif")  
    if os.path.exists(ica_filename) is False:
        ica = ICA(n_components=variance, method='fastica', max_iter=500) # , max_iter=500
        ica.fit(raw, picks=select_sensors, reject=reject) # decim = 3, 

        has_ICA = False
    else:
        has_ICA = True
        print ica_filename + '   exists!!!'
        ica = read_ica(ica_filename)
        ica.exclude = []

    # 2) identify bad components by analyzing latent sources.
    # generate ECG epochs use detection via phase statistics

    # if we just have exclude channels we jump these steps
#    if len(ica.exclude)==0:
    n_max_ecg = 3
    n_max_eog = 2

    # check if ECG_ch_name is in the raw channels
    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 default
    ecg_inds, scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
    print scores
    print '\n len ecg_inds *** ' + str(len(ecg_inds)) + '***\n'
    if len(ecg_inds) > 0:
        ecg_evoked = ecg_epochs.average()

        fig1 = ica.plot_scores(scores, exclude=ecg_inds,
#.........这里部分代码省略.........

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)

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

        print('Pre-selected comps: '+str(icacomps.exclude))
        print('##################')
        icacomps.excludeold=icacomps.exclude
        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 = '> '

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

def test_find_ecg():
    """Test find ECG peaks."""
    # Test if ECG analysis will work on data that is not preloaded
    raw = read_raw_fif(raw_fname, preload=False)

    # once with mag-trick
    # once with characteristic channel
    raw_bad = raw.copy().load_data()
    ecg_idx = raw.ch_names.index('MEG 1531')
    raw_bad._data[ecg_idx, :1] = 1e6  # this will break the detector
    raw_bad.annotations.append(raw.first_samp / raw.info['sfreq'],
                               1. / raw.info['sfreq'], 'BAD_values')
    for ch_name in ['MEG 1531', None]:
        events, ch_ECG, average_pulse, ecg = find_ecg_events(
            raw, event_id=999, ch_name=ch_name, return_ecg=True)
        assert raw.n_times == ecg.shape[-1]
        n_events = len(events)
        _, times = raw[0, :]
        assert 55 < average_pulse < 60
        # with annotations
        with pytest.deprecated_call():
            average_pulse = find_ecg_events(raw_bad, ch_name=ch_name)[2]
        assert average_pulse < 1.
        average_pulse = find_ecg_events(raw_bad, ch_name=ch_name,
                                        reject_by_annotation=True)[2]
        assert 55 < average_pulse < 60
    average_pulse = find_ecg_events(raw_bad, ch_name='MEG 2641',
                                    reject_by_annotation=False)[2]
    assert 55 < average_pulse < 65
    del raw_bad

    picks = pick_types(
        raw.info, meg='grad', eeg=False, stim=False,
        eog=False, ecg=True, emg=False, ref_meg=False,
        exclude='bads')

    # There should be no ECG channels, or else preloading will not be
    # tested
    assert 'ecg' not in raw

    ecg_epochs = create_ecg_epochs(raw, picks=picks, keep_ecg=True)
    assert len(ecg_epochs.events) == n_events
    assert 'ECG-SYN' not in raw.ch_names
    assert 'ECG-SYN' in ecg_epochs.ch_names

    picks = pick_types(
        ecg_epochs.info, meg=False, eeg=False, stim=False,
        eog=False, ecg=True, emg=False, ref_meg=False,
        exclude='bads')
    assert len(picks) == 1

    ecg_epochs = create_ecg_epochs(raw, ch_name='MEG 2641')
    assert 'MEG 2641' in ecg_epochs.ch_names

    # test with user provided ecg channel
    raw.info['projs'] = list()
    with pytest.warns(RuntimeWarning, match='unit for channel'):
        raw.set_channel_types({'MEG 2641': 'ecg'})
    create_ecg_epochs(raw)

    raw.load_data().pick_types()  # remove ECG
    ecg_epochs = create_ecg_epochs(raw, keep_ecg=False)
    assert len(ecg_epochs.events) == n_events
    assert 'ECG-SYN' not in raw.ch_names
    assert 'ECG-SYN' not in ecg_epochs.ch_names

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

# maximum number of components to reject
n_max_ecg, n_max_eog = 3, 1  # here we don't expect horizontal EOG components

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

title = 'Sources related to %s artifacts (red)'

# generate ECG epochs use detection via phase statistics

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)

def run_epochs(subject_id):
    subject = "sub%03d" % subject_id
    print("processing subject: %s" % subject)

    data_path = op.join(meg_dir, subject)

    all_epochs = list()

    # Get all bad channels
    mapping = map_subjects[subject_id]  # map to correct subject
    all_bads = list()
    for run in range(1, 7):
        bads = list()
        bad_name = op.join('bads', mapping, 'run_%02d_raw_tr.fif_bad' % run)
        if os.path.exists(bad_name):
            with open(bad_name) as f:
                for line in f:
                    bads.append(line.strip())
        all_bads += [bad for bad in bads if bad not in all_bads]

    for run in range(1, 7):
        print " - Run %s" % run
        run_fname = op.join(data_path, 'run_%02d_filt_sss_raw.fif' % run)
        if not os.path.exists(run_fname):
            continue

        raw = mne.io.Raw(run_fname, preload=True, add_eeg_ref=False)

        raw.set_channel_types({'EEG061': 'eog',
                               'EEG062': 'eog',
                               'EEG063': 'ecg',
                               'EEG064': 'misc'})  # EEG064 free floating el.
        raw.rename_channels({'EEG061': 'EOG061',
                             'EEG062': 'EOG062',
                             'EEG063': 'ECG063'})

        eog_events = mne.preprocessing.find_eog_events(raw)
        eog_events[:, 0] -= int(0.25 * raw.info['sfreq'])
        annotations = mne.Annotations(eog_events[:, 0] / raw.info['sfreq'],
                                      np.repeat(0.5, len(eog_events)),
                                      'BAD_blink', raw.info['meas_date'])
        raw.annotations = annotations  # Remove epochs with blinks

        delay = int(0.0345 * raw.info['sfreq'])
        events = mne.read_events(op.join(data_path,
                                         'run_%02d_filt_sss-eve.fif' % run))

        events[:, 0] = events[:, 0] + delay

        raw.info['bads'] = all_bads
        raw.interpolate_bads()
        raw.set_eeg_reference()

        picks = mne.pick_types(raw.info, meg=True, eeg=True, stim=True,
                               eog=True)

        # Read epochs
        epochs = mne.Epochs(raw, events, events_id, tmin, tmax, proj=True,
                            picks=picks, baseline=baseline, preload=True,
                            decim=2, reject=reject, add_eeg_ref=False)

        # ICA
        ica_name = op.join(meg_dir, subject, 'run_%02d-ica.fif' % run)
        ica = read_ica(ica_name)
        n_max_ecg = 3  # use max 3 components
        ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5)
        ecg_inds, scores_ecg = ica.find_bads_ecg(ecg_epochs, method='ctps',
                                                 threshold=0.8)
        ica.exclude += ecg_inds[:n_max_ecg]

        ica.apply(epochs)
        all_epochs.append(epochs)

    epochs = mne.epochs.concatenate_epochs(all_epochs)
    epochs.save(op.join(data_path, '%s-epo.fif' % subject))

    # To save an ICA solution you can say:

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

title = 'Sources related to %s artifacts (red)'

# generate ECG epochs use detection via phase statistics

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

# create_ecg_epochs is strange: it strips the channels of anything non M/EEG
# UNLESS picks=None
picks=None
ecg_epochs = create_ecg_epochs(raw, ch_name='ECG002', tmin=-.5, tmax=.5, picks=picks, verbose=True)

# This will work with the above, but uses MASSIV