Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Project CSP patterns to source #950

Draft
wants to merge 9 commits into
base: main
Choose a base branch
from
Draft

Project CSP patterns to source #950

Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
d15b9a8
save patterns and filters
SophieHerbst May 15, 2024
4362979
save patterns and filters - missing import
SophieHerbst May 15, 2024
6cc2e23
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] May 15, 2024
8732634
save patterns and filters - add linear model
SophieHerbst May 15, 2024
c161189
Merge branch 'main' into CSP2source
SophieHerbst May 15, 2024
814239d
complete merge
SophieHerbst May 15, 2024
a0b7a06
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] May 15, 2024
92cfef9
write csp output to proper directory
SophieHerbst May 15, 2024
a01544b
merge error
SophieHerbst May 15, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 43 additions & 5 deletions mne_bids_pipeline/steps/sensor/_05_decoding_csp.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
import mne
import numpy as np
import pandas as pd
from mne.decoding import CSP
from mne.decoding import CSP, LinearModel
from mne_bids import BIDSPath
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.pipeline import make_pipeline
Expand Down Expand Up @@ -162,6 +162,10 @@ def one_subject_decoding(
epochs=epochs,
)

# Create output directory if it doesn't already exist
output_dir = bids_path.fpath.parent / "CSD_output"
output_dir.mkdir(exist_ok=True)

# Classifier
csp = CSP(
n_components=4, # XXX revisit
Expand All @@ -170,10 +174,12 @@ def one_subject_decoding(
clf = make_pipeline(
*preproc_steps,
csp,
LogReg(
solver="liblinear", # much faster than the default
random_state=cfg.random_state,
n_jobs=1,
LinearModel(
LogReg(
solver="liblinear", # much faster than the default
random_state=cfg.random_state,
n_jobs=1,
)
),
)
cv = StratifiedKFold(
Expand Down Expand Up @@ -239,6 +245,7 @@ def _fmt_contrast(cond1, cond2, fmin, fmax, freq_range_name, tmin=None, tmax=Non
epochs_filt, y = prepare_epochs_and_y(
epochs=epochs, contrast=contrast, fmin=fmin, fmax=fmax, cfg=cfg
)

# Get the data for all time points
X = epochs_filt.get_data()

Expand All @@ -253,6 +260,21 @@ def _fmt_contrast(cond1, cond2, fmin, fmax, freq_range_name, tmin=None, tmax=Non
freq_decoding_table.loc[idx, "mean_crossval_score"] = cv_scores.mean()
freq_decoding_table.at[idx, "scores"] = cv_scores

# COEFS
clf.fit(X, y)
weights_csp = mne.decoding.get_coef(clf, "patterns_", inverse_transform=True)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Typically this is all I use. Why also take the csp.patterns_? csp.fit(X, y) results anywhere else so it's a bit weird to me to add them here. It seems like we should just use these weights_csp. And maybe we should call them clf_patterns_ because they're really the patterns inverse transformed from the CSP all the way back through the other steps (e.g., PCA, sensor scaling)?

Copy link
Collaborator Author

@SophieHerbst SophieHerbst May 15, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

that is exactly where I am lost in the terminology..

my approach was to save the weights in this step to pick them up in a separate source projection step later

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wanted to save both, because in @britta-wstnr's code they are combined during beamforming:

stc_csp = beamform_components(weights_csp, sensor_pattern_csp, spat_filter,
                                fwd, multipliers=multiplier)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The idea here is to combine the coefficients and the CSP components. The (Haufe et al.) patterns are not easily beamform-able as they do not represent classical sink-source patterns anymore but refer to the decoding classes. Using the CSP components is a work-around for that.


# PATTERNS
csp.fit_transform(X, y)
sensor_pattern_csp = csp.patterns_

# save weights and patterns
csp_patterns_fname = f"{cond1}_{cond2}_{str(fmin)}_{str(fmax)}_Hz_patterns"
csp_weights_fname = f"{cond1}_{cond2}_{str(fmin)}_{str(fmax)}_Hz_weights"

np.save(op.join(output_dir, csp_patterns_fname), sensor_pattern_csp)
np.save(op.join(output_dir, csp_weights_fname), weights_csp)

# Loop over times x frequencies
#
# Note: We don't support varying time ranges for different frequency
Expand Down Expand Up @@ -306,6 +328,7 @@ def _fmt_contrast(cond1, cond2, fmin, fmax, freq_range_name, tmin=None, tmax=Non
# Crop data to the time window of interest
if tmax is not None: # avoid warnings about outside the interval
tmax = min(tmax, epochs_filt.times[-1])

X = epochs_filt.crop(tmin, tmax).get_data()
del epochs_filt
cv_scores = cross_val_score(
Expand All @@ -323,6 +346,21 @@ def _fmt_contrast(cond1, cond2, fmin, fmax, freq_range_name, tmin=None, tmax=Non
msg += f": {cfg.decoding_metric}={score:0.3f}"
logger.info(**gen_log_kwargs(msg))

# COEFS
clf.fit(X, y)
weights_csp = mne.decoding.get_coef(clf, "patterns_", inverse_transform=True)

# PATTERNS
csp.fit_transform(X, y)
sensor_pattern_csp = csp.patterns_

# save weights and patterns
csp_patterns_fname = f"{cond1}_{cond2}_{str(fmin)}_{str(fmax)}_Hz_{str(tmin)}_{str(tmax)}_s_patterns"
csp_weights_fname = f"{cond1}_{cond2}_{str(fmin)}_{str(fmax)}_Hz_{str(tmin)}_{str(tmax)}_s_patterns"

np.save(op.join(output_dir, csp_patterns_fname), sensor_pattern_csp)
np.save(op.join(output_dir, csp_weights_fname), weights_csp)

# Write each DataFrame to a different Excel worksheet.
a_vs_b = f"{condition1}+{condition2}".replace(op.sep, "")
processing = f"{a_vs_b}+CSP+{cfg.decoding_metric}"
Expand Down
Loading