1"""
2Searchlight analysis of face vs house recognition
3==================================================
4
5Searchlight analysis requires fitting a classifier a large amount of
6times. As a result, it is an intrinsically slow method. In order to speed
7up computing, in this example, Searchlight is run only on one slice on
8the fMRI (see the generated figures).
9
10"""
11
12#########################################################################
13# Load Haxby dataset
14# -------------------
15import pandas as pd
16from nilearn import datasets
17from nilearn.image import new_img_like, load_img, get_data
18
19# We fetch 2nd subject from haxby datasets (which is default)
20haxby_dataset = datasets.fetch_haxby()
21
22# print basic information on the dataset
23print('Anatomical nifti image (3D) is located at: %s' % haxby_dataset.mask)
24print('Functional nifti image (4D) is located at: %s' % haxby_dataset.func[0])
25
26fmri_filename = haxby_dataset.func[0]
27labels = pd.read_csv(haxby_dataset.session_target[0], sep=" ")
28y = labels['labels']
29session = labels['chunks']
30
31#########################################################################
32# Restrict to faces and houses
33# ------------------------------
34from nilearn.image import index_img
35condition_mask = y.isin(['face', 'house'])
36
37fmri_img = index_img(fmri_filename, condition_mask)
38y, session = y[condition_mask], session[condition_mask]
39
40#########################################################################
41# Prepare masks
42# --------------
43# - mask_img is the original mask
44# - process_mask_img is a subset of mask_img, it contains the voxels that
45#   should be processed (we only keep the slice z = 26 and the back of the
46#   brain to speed up computation)
47import numpy as np
48
49mask_img = load_img(haxby_dataset.mask)
50
51# .astype() makes a copy.
52process_mask = get_data(mask_img).astype(np.int)
53picked_slice = 29
54process_mask[..., (picked_slice + 1):] = 0
55process_mask[..., :picked_slice] = 0
56process_mask[:, 30:] = 0
57process_mask_img = new_img_like(mask_img, process_mask)
58
59#########################################################################
60# Searchlight computation
61# -------------------------
62
63# Make processing parallel
64# /!\ As each thread will print its progress, n_jobs > 1 could mess up the
65#     information output.
66n_jobs = 1
67
68# Define the cross-validation scheme used for validation.
69# Here we use a KFold cross-validation on the session, which corresponds to
70# splitting the samples in 4 folds and make 4 runs using each fold as a test
71# set once and the others as learning sets
72from sklearn.model_selection import KFold
73cv = KFold(n_splits=4)
74
75import nilearn.decoding
76# The radius is the one of the Searchlight sphere that will scan the volume
77searchlight = nilearn.decoding.SearchLight(
78    mask_img,
79    process_mask_img=process_mask_img,
80    radius=5.6, n_jobs=n_jobs,
81    verbose=1, cv=cv)
82searchlight.fit(fmri_img, y)
83
84#########################################################################
85# F-scores computation
86# ----------------------
87from nilearn.input_data import NiftiMasker
88
89# For decoding, standardizing is often very important
90nifti_masker = NiftiMasker(mask_img=mask_img, sessions=session,
91                           standardize=True, memory='nilearn_cache',
92                           memory_level=1)
93fmri_masked = nifti_masker.fit_transform(fmri_img)
94
95from sklearn.feature_selection import f_classif
96f_values, p_values = f_classif(fmri_masked, y)
97p_values = -np.log10(p_values)
98p_values[p_values > 10] = 10
99p_unmasked = get_data(nifti_masker.inverse_transform(p_values))
100
101#########################################################################
102# Visualization
103# --------------
104# Use the fmri mean image as a surrogate of anatomical data
105from nilearn import image
106mean_fmri = image.mean_img(fmri_img)
107
108from nilearn.plotting import plot_stat_map, plot_img, show
109searchlight_img = new_img_like(mean_fmri, searchlight.scores_)
110
111# Because scores are not a zero-center test statistics, we cannot use
112# plot_stat_map
113plot_img(searchlight_img, bg_img=mean_fmri,
114         title="Searchlight", display_mode="z", cut_coords=[-9],
115         vmin=.42, cmap='hot', threshold=.2, black_bg=True)
116
117# F_score results
118p_ma = np.ma.array(p_unmasked, mask=np.logical_not(process_mask))
119f_score_img = new_img_like(mean_fmri, p_ma)
120plot_stat_map(f_score_img, mean_fmri,
121              title="F-scores", display_mode="z",
122              cut_coords=[-9],
123              colorbar=False)
124
125show()
126