Python decoding.GeneralizationAcrossTime类代码示例

def quick_score(X, y, clf=None, scorer=None):
    from sklearn.cross_validation import KFold
    regression = (len(np.unique(y)) > 2) & isinstance(y[0], float)
    if scorer is None:
        scorer = scorer_spearman if regression else scorer_auc
    if clf is None:
        clf = RidgeCV(alphas=[(2 * C) ** -1 for C in [1e-4, 1e-2, 1]])\
            if regression else force_predict(LogisticRegression(), axis=1)
    sel = np.where(~np.isnan(y))[0]
    X = X[sel, :, :]
    y = y[sel]
    epochs = mat2mne(X, sfreq=100)
    clf = make_pipeline(StandardScaler(), clf)
    cv = KFold(len(y), 5) if regression else None
    gat = GeneralizationAcrossTime(clf=clf, n_jobs=-1, scorer=scorer, cv=cv)
    gat.fit(epochs, y)
    gat.score(epochs, y)
    return gat

def _get_data(tmin=-0.2, tmax=0.5, event_id=dict(aud_l=1, vis_l=3),
              event_id_gen=dict(aud_l=2, vis_l=4), test_times=None):
    """Aux function for testing GAT viz."""
    with warnings.catch_warnings(record=True):  # deprecated
        gat = GeneralizationAcrossTime()
    raw = read_raw_fif(raw_fname)
    raw.add_proj([], remove_existing=True)
    events = read_events(event_name)
    picks = pick_types(raw.info, meg='mag', stim=False, ecg=False,
                       eog=False, exclude='bads')
    picks = picks[1:13:3]
    decim = 30
    # Test on time generalization within one condition
    with warnings.catch_warnings(record=True):
        epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                        preload=True, decim=decim)
    epochs_list = [epochs[k] for k in event_id]
    equalize_epoch_counts(epochs_list)
    epochs = concatenate_epochs(epochs_list)

    # Test default running
    with warnings.catch_warnings(record=True):  # deprecated
        gat = GeneralizationAcrossTime(test_times=test_times)
    gat.fit(epochs)
    gat.score(epochs)
    return gat

def _run(epochs, events, analysis):
    """Runs temporal generalization for a given subject and analysis"""
    print(subject, analysis['name'])

    # subselect the trials (e.g. exclude absent trials) with a
    # dataframe query defined in conditions.py
    query, condition = analysis['query'], analysis['condition']
    sel = range(len(events)) if query is None \
        else events.query(query).index
    sel = [ii for ii in sel if ~np.isnan(events[condition][sel][ii])]

    # The to-be-predicted value, for each trial:
    y = np.array(events[condition], dtype=np.float32)

    print analysis['name'], np.unique(y[sel]), len(sel)

    # Abort if there is no trial
    if len(sel) == 0:
        return

    # Apply analysis
    gat = GeneralizationAcrossTime(clf=analysis['clf'],
                                   cv=analysis['cv'],
                                   scorer=analysis['scorer'],
                                   n_jobs=-1)
    print(subject, analysis['name'], 'fit')
    gat.fit(epochs[sel], y=y[sel])
    print(subject, analysis['name'], 'score')
    score = gat.score(epochs[sel], y=y[sel])
    print(subject, analysis['name'], 'save')

    # save space
    if analysis['name'] not in ['probe_phase', 'target_circAngle']:
        # we'll need the estimator trained on the probe_phase and to generalize
        # to the target phase and prove that there is a significant signal.
        gat.estimators_ = None
    if analysis['name'] not in ['target_present', 'target_circAngle',
                                'probe_circAngle']:
        # We need these individual prediction to control for the correlation
        # between target and probe angle.
        gat.y_pred_ = None

    # Save analysis
    save([gat, analysis, sel, events], 'decod',
         subject=subject, analysis=analysis['name'], overwrite=True,
         upload=True)
    save([score, epochs.times], 'score',
         subject=subject, analysis=analysis['name'], overwrite=True,
         upload=True)
    return

def _decod(subject, analysis):
    from mne.decoding import GeneralizationAcrossTime

    # if already computed let's just load it from disk
    fname_kwargs = dict(subject=subject, analysis=analysis['name'] + '_vhp')
    score_fname = paths('score', **fname_kwargs)
    if op.exists(score_fname):
        return load('score', **fname_kwargs)

    epochs = _get_epochs(subject)
    events = load('behavior', subject=subject)

    # Let's not recompute everything, this is just a control analysis
    print(subject, analysis['name'])
    epochs._data = epochs.get_data()
    epochs.preload = True
    epochs.crop(0., .900)
    epochs.decimate(2)

    query, condition = analysis['query'], analysis['condition']
    sel = range(len(events)) if query is None else events.query(query).index
    sel = [ii for ii in sel if ~np.isnan(events[condition][sel][ii])]
    y = np.array(events[condition], dtype=np.float32)

    print analysis['name'], np.unique(y[sel]), len(sel)

    if len(sel) == 0:
        return

    # Apply analysis
    gat = GeneralizationAcrossTime(clf=analysis['clf'],
                                   cv=analysis['cv'],
                                   scorer=analysis['scorer'],
                                   n_jobs=-1)
    print(subject, analysis['name'], 'fit')
    gat.fit(epochs[sel], y=y[sel])
    print(subject, analysis['name'], 'score')
    score = gat.score(epochs[sel], y=y[sel])
    print(subject, analysis['name'], 'save')

    # save space
    gat.estimators_ = None
    gat.y_pred_ = None

    # Save analysis
    save([score, epochs.times], 'score', overwrite=True, upload=True,
         **fname_kwargs)
    return score, epochs.times

def _get_data():
    """Aux function for testing GAT viz"""
    gat = GeneralizationAcrossTime()
    raw = io.Raw(raw_fname, preload=False)
    events = read_events(event_name)
    picks = pick_types(raw.info, meg='mag', stim=False, ecg=False,
                       eog=False, exclude='bads')
    picks = picks[1:13:3]
    decim = 30
    # Test on time generalization within one condition
    with warnings.catch_warnings(record=True):
        epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                        baseline=(None, 0), preload=True, decim=decim)

    # Test default running
    gat = GeneralizationAcrossTime()
    gat.fit(epochs)
    gat.score(epochs)
    return gat

info = create_info(chan_names, 1, chan_types)
events = np.c_[np.cumsum(np.ones(n_trial)), np.zeros(n_trial), np.zeros(n_trial)]
epochs = EpochsArray(data, info, events)

# RUN GAT ======================================================================

# SVR
# --- fit & predict separately
cos = lambda angles: np.cos(angle2circle(angles))
sin = lambda angles: np.sin(angle2circle(angles))
gats = list()
for transform in [cos, sin]:
    scaler = StandardScaler()
    svr = SVR(C=1, kernel='linear')
    clf = Pipeline([('scaler', scaler), ('svr', svr)])
    gat = GeneralizationAcrossTime(n_jobs=-1, clf=clf)
    gat.fit(epochs, y=transform(trial_angles))
    gat.predict(epochs)
    gats.append(gat)
# --- recombine
predict_angles, true_angles = recombine_svr_prediction(gats[0], gats[1])
# --- score
angle_errors_svr = compute_error_svr(predict_angles, true_angles)
plt.matshow(np.mean(angle_errors_svr,axis=2)), plt.colorbar(), plt.show()


# SVC Gat
scaler = StandardScaler()
svc = SVC(C=1, kernel='linear', probability=True)
clf = Pipeline([('scaler', scaler), ('svc', svc)])
gat = GeneralizationAcrossTime(n_jobs=-1, clf=clf, predict_type='predict_proba')

# Create epochs to use for classification
n_trial, n_chan, n_time = X.shape
events = np.vstack((range(n_trial), np.zeros(n_trial, int), y.astype(int))).T
chan_names = ['MEG %i' % chan for chan in range(n_chan)]
chan_types = ['mag'] * n_chan
sfreq = 250
info = create_info(chan_names, sfreq, chan_types)
epochs = EpochsArray(data=X, info=info, events=events, verbose=False)
epochs.times = selected_times[:n_time]

# make classifier
clf = LogisticRegression(C=0.0001)

# fit model and score
gat = GeneralizationAcrossTime(
    clf=clf, scorer="roc_auc", cv=cv, predict_method="predict")
gat.fit(epochs, y=y)
gat.score(epochs, y=y)

# Save model
joblib.dump(gat, data_path + "decode_time_gen/gat_ge.jl")

# make matrix plot and save it
fig = gat.plot(
    cmap="viridis", title="Temporal Gen (Classic vs planning) for Global Eff.")
fig.savefig(data_path + "decode_time_gen/gat_matrix_ge.png")

fig = gat.plot_diagonal(
    chance=0.5, title="Temporal Gen (Classic vs planning) for Global eff.")
fig.savefig(data_path + "decode_time_gen/gat_diagonal_ge.png")

def test_circular_classifiers():
    from mne.decoding import GeneralizationAcrossTime
    from ..scorers import scorer_angle
    from sklearn.linear_model import Ridge, RidgeCV
    epochs, angles = make_circular_data()
    clf_list = [PolarRegression, AngularRegression,
                SVR_polar, SVR_angle]  # XXX will be deprecated
    for clf_init in clf_list:
        for independent in [False, True]:
            if clf_init in [SVR_polar, SVR_angle]:
                if (not independent):
                    continue
                clf = clf_init(clf=Ridge(random_state=0))
            else:
                clf = clf_init(clf=Ridge(random_state=0),
                               independent=independent)
            print clf_init, independent
            gat = GeneralizationAcrossTime(clf=clf, scorer=scorer_angle)
            gat.fit(epochs, y=angles)
            gat.predict(epochs)
            gat.score(y=angles)
            assert_true(np.abs(gat.scores_[0][0]) < .5)  # chance level
            assert_true(gat.scores_[1][1] > 1.)  # decode
            assert_true(gat.scores_[2][2] > 1.)  # decode
            assert_true(gat.scores_[1][2] < -1.)  # anti-generalize
    # Test args
    gat = GeneralizationAcrossTime(clf=RidgeCV(alphas=[1., 2.]),
                                   scorer=scorer_angle)
    gat.fit(epochs, y=angles)
    gat = GeneralizationAcrossTime(clf=RidgeCV(), scorer=scorer_angle)
    gat.fit(epochs, y=angles)

cv = StratifiedKFold(n_splits=10, shuffle=True)

# Create epochs to use for classification
n_trial, n_chan, n_time = X.shape
events = np.vstack((range(n_trial), np.zeros(n_trial, int), y.astype(int))).T
chan_names = ['MEG %i' % chan for chan in range(n_chan)]
chan_types = ['mag'] * n_chan
sfreq = 250
info = create_info(chan_names, sfreq, chan_types)
epochs = EpochsArray(data=X, info=info, events=events, verbose=False)
epochs.times = selected_times[:n_time]

epochs.crop(-3.8, None)

# fit model and score
gat = GeneralizationAcrossTime(
    scorer="accuracy", cv=cv, predict_method="predict")
gat.fit(epochs, y=y)
gat.score(epochs, y=y)

# Save model
joblib.dump(gat, data_path + "decode_time_gen/%s_gat_tr.jl" % subject)

# make matrix plot and save it
fig = gat.plot(
    cmap="viridis",
    title="Temporal Gen (Classic vs planning) for transitivity.")
fig.savefig(data_path + "decode_time_gen/%s_gat_matrix_tr.png" % subject)

fig = gat.plot_diagonal(
    chance=0.5, title="Temporal Gen (Classic vs planning) for transitivity")
fig.savefig(data_path + "decode_time_gen/%s_gat_diagonal_tr.png" % subject)

def test_generalization_across_time():
    """Test time generalization decoding
    """
    from sklearn.svm import SVC
    from sklearn.preprocessing import LabelEncoder
    from sklearn.metrics import mean_squared_error

    raw = io.Raw(raw_fname, preload=False)
    events = read_events(event_name)
    picks = pick_types(raw.info, meg='mag', stim=False, ecg=False,
                       eog=False, exclude='bads')
    picks = picks[0:2]
    decim = 30

    # Test on time generalization within one condition
    with warnings.catch_warnings(record=True):
        epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks,
                        baseline=(None, 0), preload=True, decim=decim)
    # Test default running
    gat = GeneralizationAcrossTime()
    assert_equal("<GAT | no fit, no prediction, no score>", "%s" % gat)
    assert_raises(ValueError, gat.fit, epochs, picks='foo')
    with warnings.catch_warnings(record=True):
        # check classic fit + check manual picks
        gat.fit(epochs, picks=[0])
        # check optional y as array
        gat.fit(epochs, y=epochs.events[:, 2])
        # check optional y as list
        gat.fit(epochs, y=epochs.events[:, 2].tolist())
    assert_equal(len(gat.picks_), len(gat.ch_names), 1)
    assert_equal("<GAT | fitted, start : -0.200 (s), stop : 0.499 (s), no "
                 "prediction, no score>", '%s' % gat)
    assert_equal(gat.ch_names, epochs.ch_names)
    gat.predict(epochs)
    assert_equal("<GAT | fitted, start : -0.200 (s), stop : 0.499 (s), "
                 "predicted 14 epochs, no score>",
                 "%s" % gat)
    gat.score(epochs)
    gat.score(epochs, y=epochs.events[:, 2])
    gat.score(epochs, y=epochs.events[:, 2].tolist())
    assert_equal("<GAT | fitted, start : -0.200 (s), stop : 0.499 (s), "
                 "predicted 14 epochs,\n scored "
                 "(accuracy_score)>", "%s" % gat)
    with warnings.catch_warnings(record=True):
        gat.fit(epochs, y=epochs.events[:, 2])

    old_mode = gat.predict_mode
    gat.predict_mode = 'super-foo-mode'
    assert_raises(ValueError, gat.predict, epochs)
    gat.predict_mode = old_mode

    gat.score(epochs, y=epochs.events[:, 2])
    assert_true("accuracy_score" in '%s' % gat.scorer_)
    epochs2 = epochs.copy()

    # check _DecodingTime class
    assert_equal("<DecodingTime | start: -0.200 (s), stop: 0.499 (s), step: "
                 "0.047 (s), length: 0.047 (s), n_time_windows: 15>",
                 "%s" % gat.train_times)
    assert_equal("<DecodingTime | start: -0.200 (s), stop: 0.499 (s), step: "
                 "0.047 (s), length: 0.047 (s), n_time_windows: 15 x 15>",
                 "%s" % gat.test_times_)

    # the y-check
    gat.predict_mode = 'mean-prediction'
    epochs2.events[:, 2] += 10
    gat_ = copy.deepcopy(gat)
    assert_raises(ValueError, gat_.score, epochs2)
    gat.predict_mode = 'cross-validation'

    # Test basics
    # --- number of trials
    assert_true(gat.y_train_.shape[0] ==
                gat.y_true_.shape[0] ==
                len(gat.y_pred_[0][0]) == 14)
    # ---  number of folds
    assert_true(np.shape(gat.estimators_)[1] == gat.cv)
    # ---  length training size
    assert_true(len(gat.train_times['slices']) == 15 ==
                np.shape(gat.estimators_)[0])
    # ---  length testing sizes
    assert_true(len(gat.test_times_['slices']) == 15 ==
                np.shape(gat.scores_)[0])
    assert_true(len(gat.test_times_['slices'][0]) == 15 ==
                np.shape(gat.scores_)[1])

    # Test longer time window
    gat = GeneralizationAcrossTime(train_times={'length': .100})
    with warnings.catch_warnings(record=True):
        gat2 = gat.fit(epochs)
    assert_true(gat is gat2)  # return self
    assert_true(hasattr(gat2, 'cv_'))
    assert_true(gat2.cv_ != gat.cv)
    scores = gat.score(epochs)
    assert_true(isinstance(scores, list))  # type check
    assert_equal(len(scores[0]), len(scores))  # shape check

    assert_equal(len(gat.test_times_['slices'][0][0]), 2)
    # Decim training steps
    gat = GeneralizationAcrossTime(train_times={'step': .100})
#.........这里部分代码省略.........

###############################################################################
# Generalization Across Time
# --------------------------
#
# This runs the analysis used in [1]_ and further detailed in [2]_
#
# Here we'll use a stratified cross-validation scheme.

# make response vector
y = np.zeros(len(epochs.events), dtype=int)
y[epochs.events[:, 2] == 3] = 1
cv = StratifiedKFold(y=y)  # do a stratified cross-validation

# define the GeneralizationAcrossTime object
gat = GeneralizationAcrossTime(predict_mode="cross-validation", n_jobs=1, cv=cv, scorer=roc_auc_score)

# fit and score
gat.fit(epochs, y=y)
gat.score(epochs)

# let's visualize now
gat.plot()
gat.plot_diagonal()

###############################################################################
# Exercise
# --------
#  - Can you improve the performance using full epochs and a common spatial
#    pattern (CSP) used by most BCI systems?
#  - Explore other datasets from MNE (e.g. Face dataset from SPM to predict

decim = 2  # decimate to make the example faster to run
epochs = mne.Epochs(raw, events, event_id, -0.050, 0.400, proj=True,
                    picks=picks, baseline=None, preload=True,
                    reject=dict(mag=5e-12), decim=decim, verbose=False)

# We will train the classifier on all left visual vs auditory trials
# and test on all right visual vs auditory trials.

# In this case, because the test data is independent from the train data,
# we test the classifier of each fold and average the respective predictions.

# Define events of interest
triggers = epochs.events[:, 2]
viz_vs_auditory = np.in1d(triggers, (1, 2)).astype(int)

gat = GeneralizationAcrossTime(predict_mode='mean-prediction', n_jobs=1)

# For our left events, which ones are visual?
viz_vs_auditory_l = (triggers[np.in1d(triggers, (1, 3))] == 3).astype(int)
# To make scikit-learn happy, we converted the bool array to integers
# in the same line. This results in an array of zeros and ones:
print("The unique classes' labels are: %s" % np.unique(viz_vs_auditory_l))

gat.fit(epochs[('AudL', 'VisL')], y=viz_vs_auditory_l)

# For our right events, which ones are visual?
viz_vs_auditory_r = (triggers[np.in1d(triggers, (2, 4))] == 4).astype(int)

gat.score(epochs[('AudR', 'VisR')], y=viz_vs_auditory_r)
gat.plot(
    title="Generalization Across Time (visual vs auditory): left to right")

            query, condition = analysis['query'], analysis['condition']
            sel = range(len(events)) if query is None \
                else events.query(query).index
            sel = [ii for ii in sel if ~np.isnan(events[condition][sel][ii])]
            y = np.array(events[condition], dtype=np.float32)

            print analysis['name'], np.unique(y[sel]), len(sel)

            if len(sel) == 0:
                logger.warning('%s: no epoch in %s for %s.' % (
                    subject, data_type, analysis['name']))
                continue

            # Apply analysis
            gat = GeneralizationAcrossTime(clf=analysis['clf'],
                                           cv=analysis['cv'],
                                           scorer=analysis['scorer'],
                                           n_jobs=-1)
            gat.fit(epochs[sel], y=y[sel])
            gat.score(epochs[sel], y=y[sel])

            # Save analysis
            pkl_fname = paths('decod', subject=subject, data_type=data_type,
                              analysis=analysis['name'], log=True)

            # Save classifier results
            with open(pkl_fname, 'wb') as f:
                pickle.dump([gat, analysis, sel, events], f)

            # Plot
            fig = gat.plot_diagonal(show=False)
            report.add_figs_to_section(fig, ('%s %s %s: (diagonal)' %

print(__doc__)

# Preprocess data
data_path = spm_face.data_path()
# Load and filter data, set up epochs
raw_fname = data_path + '/MEG/spm/SPM_CTF_MEG_example_faces%d_3D_raw.fif'

raw = mne.io.Raw(raw_fname % 1, preload=True)  # Take first run

picks = mne.pick_types(raw.info, meg=True, exclude='bads')
raw.filter(1, 45, method='iir')

events = mne.find_events(raw, stim_channel='UPPT001')
event_id = {"faces": 1, "scrambled": 2}
tmin, tmax = -0.1, 0.5

decim = 4  # decimate to make the example faster to run
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
                    picks=picks, baseline=None, preload=True,
                    reject=dict(mag=1.5e-12), decim=decim, verbose=False)

# Define decoder. The decision function is employed to use cross-validation
gat = GeneralizationAcrossTime(predict_mode='cross-validation', n_jobs=1)

# fit and score
gat.fit(epochs)
gat.score(epochs)
gat.plot(vmin=0.1, vmax=0.9,
         title="Generalization Across Time (faces vs. scrambled)")
gat.plot_diagonal()  # plot decoding across time (correspond to GAT diagonal)

	# decoding
	epochs = epochs.decimate(decim)
	n_tr = len(epochs.times)
	X = epochs[stim_list>0]
	y = stim_list[stim_list>0]

	# SVM parameters
	scaler = StandardScaler() # centers data by removing the mean and scales to unit variance
	# model = svm.SVC(C=1, kernel='linear', class_weight='auto')
	model = svm.LinearSVC(C=1, multi_class='ovr', class_weight='auto')
	clf = make_pipeline(scaler, model)
	cv = StratifiedKFold(y, n_fold)
	gat = GeneralizationAcrossTime(
		cv=cv,
		clf=clf, 
		predict_mode='cross-validation', 
		n_jobs=n_jobs,
		# train_times=dict(step=step,	length=length)
		)

	gat.fit(X, y)
	# score = gat.score(X, y)
	# gat.plot(vmin=.2, vmax=.3)
	# gat.plot_diagonal(chance=.25)

	#######################################################################################################
	 # now for the rsvp data 
	fname = op.join(data_path, 'abse_' + subject + '_main.mat')

	epochs = sm_fieldtrip2mne(fname)
	n_trial = len(epochs)