"""
Effect of standardization on decoding performance
=================================================

Different voxels have different baselines and variances, so the
``standardize`` parameter of nilearn maskers rescales each voxel's time
series before it reaches the estimator:

- ``"zscore_sample"`` shifts the signal to zero mean and scales it to
  unit variance,
- ``"psc"`` scales it to percent signal change relative to its mean,
- ``None`` leaves the signal as acquired.

This choice matters because estimators such as SVMs and other regularized
linear models are sensitive to feature scale. The scikit-learn section on
`standardization <https://scikit-learn.org/stable/modules/preprocessing.html#standardization-or-standard-scaling>`_
explains why a high-variance feature can dominate the objective function.

Note that the defaults differ across estimators: ``NiftiMasker`` leaves
the data untouched (``standardize=None``), while ``Decoder`` standardizes
by default (``standardize="zscore_sample"``). The code below runs the
same decoding pipeline with each setting and compares the scores.

Aside from ``standardize``, the data processing and ``Decoder``
options mirror those of the
:ref:`sphx_glr_auto_examples_02_decoding_plot_haxby_different_estimators.py`
example, so the two can be read side by side.
"""

# %%
# Load the Haxby data
# -------------------
import numpy as np
import pandas as pd

from nilearn import datasets
from nilearn.image import index_img

# by default 2nd subject data will be fetched
haxby_dataset = datasets.fetch_haxby()

fmri_filename = haxby_dataset.func[0]
mask_filename = haxby_dataset.mask_vt[0]

labels = pd.read_csv(haxby_dataset.session_target[0], sep=" ")
stimuli = labels["labels"]

# remove resting-state runs from the data and the target
task_mask = stimuli != "rest"
run_labels = labels["chunks"][task_mask]

fmri_niimgs = index_img(fmri_filename, task_mask)
classification_target = stimuli[task_mask]

# %%
# Why do the defaults differ?
# ---------------------------
# ``NiftiMasker`` is a low-level building block: it only extracts
# voxels from the mask and leaves any rescaling decision to the user,
# so it defaults to ``standardize=None``. ``Decoder`` is an
# end-to-end estimator whose regularized classifier benefits from
# rescaled features, so it applies ``standardize="zscore_sample"``
# out of the box. Rather than hard-coding them, the code below reads
# both defaults from the signatures:
import inspect

from nilearn.decoding import Decoder
from nilearn.maskers import NiftiMasker

print(
    "NiftiMasker default:",
    inspect.signature(NiftiMasker).parameters["standardize"].default,
)
print(
    "Decoder default:",
    inspect.signature(Decoder).parameters["standardize"].default,
)

# %%
# Decoding with each standardization setting
# ------------------------------------------
from sklearn.model_selection import LeaveOneGroupOut

cv = LeaveOneGroupOut()

scores = {}
for standardize in (None, "zscore_sample", "psc"):
    # everything but "standardize" is held fixed
    decoder = Decoder(
        estimator="svc",
        mask=mask_filename,
        cv=cv,
        standardize=standardize,
        screening_percentile=100,
        n_jobs=-1,
        verbose=1,
    )
    decoder.fit(
        fmri_niimgs,
        classification_target,
        groups=run_labels,
    )
    # Decoder stores one cross-validated roc_auc score per category
    scores[standardize] = np.mean(
        list(decoder.cv_scores_.values()), axis=0
    ).mean()

# %%
# Compare the scores
# ------------------
import matplotlib.pyplot as plt

from nilearn.plotting import show

setting_labels = {
    None: "standardize=None",
    "zscore_sample": 'standardize="zscore_sample"',
    "psc": 'standardize="psc"',
}

plt.figure(figsize=(7, 4))
plt.bar(
    [setting_labels[s] for s in scores],
    list(scores.values()),
    color=["C0", "C1", "C2"],
)
plt.ylim(0.5, 1.0)
plt.ylabel("cross-validated ROC AUC")
plt.title("Effect of standardization on decoding accuracy")
plt.xticks(rotation=10, ha="right")
plt.tight_layout()

show()

# %%
# What does this tell us?
# -----------------------
# The setting has a real effect: on this dataset the unstandardized
# pipeline scores highest, followed by ``"zscore_sample"`` and then
# ``"psc"``.
# This is one dataset and one classifier, so the ranking should not
# be generalized to other data. What does generalize is that the
# choice matters — and that ``"zscore_sample"`` remains a good
# default even though it does not translate into higher accuracy
# here: it puts every voxel on a common zero-mean, unit-variance
# scale, so no voxel can dominate the estimator purely because of
# its raw amplitude.

# sphinx_gallery_dummy_images=1