Nilearn: Decoder object failing weirdly when screening percentile too low.

[thomasbazeille]

On Nilearn master, running the Decoder object on Haxby example (see code below) encounters a weird bug. With screening_percentile=20, it runs smoothly but with screening_percentile=10 or lower it yields a completely unrelated error.

I don't get what is the bug here. Any ideas @tbng ?

from nilearn.datasets import fetch_haxby
data_files = fetch_haxby()

# Load behavioral data
import pandas as pd
behavioral = pd.read_csv(data_files.session_target[0], sep=" ")

# Restrict to face and house conditions
conditions = behavioral['labels']
condition_mask = conditions.isin(['face', 'house'])

# Split data into train and test samples, using the chunks
condition_mask_train = (condition_mask) & (behavioral['chunks'] <= 6)
condition_mask_test = (condition_mask) & (behavioral['chunks'] > 6)

# Apply this sample mask to X (fMRI data) and y (behavioral labels)
# Because the data is in one single large 4D image, we need to use
# index_img to do the split easily
from nilearn.image import index_img
func_filenames = data_files.func[0]
X_train = index_img(func_filenames, condition_mask_train)
X_test = index_img(func_filenames, condition_mask_test)
y_train = conditions[condition_mask_train].values
y_test = conditions[condition_mask_test].values


# Compute the mean epi to be used for the background of the plotting
from nilearn.image import mean_img
background_img = mean_img(func_filenames)

##############################################################################
# Fit fREM
# --------------------------------------
from nilearn.decoding import fREMClassifier
from nilearn.decoding import Decoder
dec = Decoder(screening_percentile=5)
dec.fit(X_train, y_train)

Full traceback :

ValueError                                Traceback (most recent call last)
<ipython-input-135-b9478a772005> in <module>()
----> 1 dec.fit(X_train, y_train)

~/Documents/Stage/nilearn/nilearn/decoding/decoder.py in fit(self, X, y, groups)
    496                 c, self.screening_percentile_, self.clustering_percentile)
    497             for c, (train, test) in itertools.product(
--> 498                 range(n_problems), self.cv_))
    499 
    500         coefs, intercepts = self._fetch_parallel_fit_outputs(

~/anaconda2/envs/py36/lib/python3.6/site-packages/joblib/parallel.py in __call__(self, iterable)
    919             # remaining jobs.
    920             self._iterating = False
--> 921             if self.dispatch_one_batch(iterator):
    922                 self._iterating = self._original_iterator is not None
    923 

~/anaconda2/envs/py36/lib/python3.6/site-packages/joblib/parallel.py in dispatch_one_batch(self, iterator)
    757                 return False
    758             else:
--> 759                 self._dispatch(tasks)
    760                 return True
    761 

~/anaconda2/envs/py36/lib/python3.6/site-packages/joblib/parallel.py in _dispatch(self, batch)
    714         with self._lock:
    715             job_idx = len(self._jobs)
--> 716             job = self._backend.apply_async(batch, callback=cb)
    717             # A job can complete so quickly than its callback is
    718             # called before we get here, causing self._jobs to

~/anaconda2/envs/py36/lib/python3.6/site-packages/joblib/_parallel_backends.py in apply_async(self, func, callback)
    180     def apply_async(self, func, callback=None):
    181         """Schedule a func to be run"""
--> 182         result = ImmediateResult(func)
    183         if callback:
    184             callback(result)

~/anaconda2/envs/py36/lib/python3.6/site-packages/joblib/_parallel_backends.py in __init__(self, batch)
    547         # Don't delay the application, to avoid keeping the input
    548         # arguments in memory
--> 549         self.results = batch()
    550 
    551     def get(self):

~/anaconda2/envs/py36/lib/python3.6/site-packages/joblib/parallel.py in __call__(self)
    223         with parallel_backend(self._backend, n_jobs=self._n_jobs):
    224             return [func(*args, **kwargs)
--> 225                     for func, args, kwargs in self.items]
    226 
    227     def __len__(self):

~/anaconda2/envs/py36/lib/python3.6/site-packages/joblib/parallel.py in <listcomp>(.0)
    223         with parallel_backend(self._backend, n_jobs=self._n_jobs):
    224             return [func(*args, **kwargs)
--> 225                     for func, args, kwargs in self.items]
    226 
    227     def __len__(self):

~/anaconda2/envs/py36/lib/python3.6/site-packages/joblib/memory.py in __call__(self, *args, **kwargs)
    353 
    354     def __call__(self, *args, **kwargs):
--> 355         return self.func(*args, **kwargs)
    356 
    357     def call_and_shelve(self, *args, **kwargs):

~/Documents/Stage/nilearn/nilearn/decoding/decoder.py in _parallel_fit(estimator, X, y, train, test, param_grid, is_classification, scorer, mask_img, class_index, screening_percentile, clustering_percentile)
    180     for param in ParameterGrid(param_grid):
    181         estimator = clone(estimator).set_params(**param)
--> 182         estimator.fit(X_train, y_train)
    183 
    184         if is_classification:

~/anaconda2/envs/py36/lib/python3.6/site-packages/sklearn/svm/classes.py in fit(self, X, y, sample_weight)
    227         X, y = check_X_y(X, y, accept_sparse='csr',
    228                          dtype=np.float64, order="C",
--> 229                          accept_large_sparse=False)
    230         check_classification_targets(y)
    231         self.classes_ = np.unique(y)

~/anaconda2/envs/py36/lib/python3.6/site-packages/sklearn/utils/validation.py in check_X_y(X, y, accept_sparse, accept_large_sparse, dtype, order, copy, force_all_finite, ensure_2d, allow_nd, multi_output, ensure_min_samples, ensure_min_features, y_numeric, warn_on_dtype, estimator)
    717                     ensure_min_features=ensure_min_features,
    718                     warn_on_dtype=warn_on_dtype,
--> 719                     estimator=estimator)
    720     if multi_output:
    721         y = check_array(y, 'csr', force_all_finite=True, ensure_2d=False,

~/anaconda2/envs/py36/lib/python3.6/site-packages/sklearn/utils/validation.py in check_array(array, accept_sparse, accept_large_sparse, dtype, order, copy, force_all_finite, ensure_2d, allow_nd, ensure_min_samples, ensure_min_features, warn_on_dtype, estimator)
    556                              " a minimum of %d is required%s."
    557                              % (n_features, array.shape, ensure_min_features,
--> 558                                 context))
    559 
    560     if warn_on_dtype and dtype_orig is not None and array.dtype != dtype_orig:

ValueError: Found array with 0 feature(s) (shape=(112, 0)) while a minimum of 1 is required.

[jeromedockes]

there seems to be a bug in the decoder: the adjust_screening_percentile is done twice, once in the decoder's fit and once in _parallel_fit, when check_feature_screening is called.
also, can you check if the masking is working well? when I run the script you posted the masking seems to fail for me: almost all the image (not just the brain) is kept inside the mask (as a result the adjustment reduces the screening percentile a lot but that's not the biggest problem)

[jeromedockes]

BTW, I think the documentation for screening_percentile could be improved a lot:

    screening_percentile: int, float, optional, in the closed interval [0, 100]
        Perform an univariate feature selection based on the Anova F-value for
        the input data. A float according to a percentile of the highest
        scores. Default: 20.

it should state that it is the percentage of brain volume that is kept, and that is expressed wrt a full mni template brain, and adjusted if the volume of your brain mask differs (for example what happens when you provide a mask selecting an ROI?)

[thomasbazeille]

Maybe related : #854

[jeromedockes]

the bug is described above. do you want to open a PR? or should I do it?

[tbng]

On 28.04.2020 06:24, jeromedockes wrote:

the bug is described above. do you want to open a PR? or should I do it?

Sorry I am a bit busy these days. I will try to investigate and open a PR tomorrow.

--
You are receiving this because you were mentioned.
Reply to this email directly or view it on GitHub:
https://github.com/nilearn/nilearn/issues/2414#issuecomment-620605798

[jeromedockes]

Sorry I am a bit busy these days. I will try to investigate and open a PR tomorrow.

no rush!

[tbng]

I can confirm that this bug can be replicated this bug with my machine with the script.

there seems to be a bug in the decoder: the adjust_screening_percentile is done twice, once in the decoder's fit and once in _parallel_fit, when check_feature_screening is called.

Good catch, but it is actually the _adjust_screening_percentile function get called twice? Once in Decoder.fit() and one more time in _parallel_fit -- a nested call for this function. This function seems to check the volume of the mask and adjust the screening percentile on that, so it is possible that the screening_percentile if inputed as very small (like 0.05 in the above script) get corrected twice might be even smaller.

also, can you check if the masking is working well? when I run the script you posted the masking seems to fail for me: almost all the image (not just the brain) is kept inside the mask (as a result the adjustment reduces the screening percentile a lot but that's not the biggest problem)

Also related that the _adjust_screening_percentile takes mask image as one of the input, so the quality of the mask has direct influence on the output.

[tbng]

BTW, I think the documentation for screening_percentile could be improved a lot:

    screening_percentile: int, float, optional, in the closed interval [0, 100]
        Perform an univariate feature selection based on the Anova F-value for
        the input data. A float according to a percentile of the highest
        scores. Default: 20.

it should state that it is the percentage of brain volume that is kept, and that is expressed wrt a full mni template brain, and adjusted if the volume of your brain mask differs (for example what happens when you provide a mask selecting an ROI?)

Yes as least we can improve this one quickly in the PR.

[thomasbazeille]

linked with #2360