From 5087a661a594af15d721d4200164a3789afc1c7b Mon Sep 17 00:00:00 2001 From: hrichard Date: Mon, 29 Nov 2021 15:25:50 +0100 Subject: [PATCH 1/2] add thomas work with fastsrm imports --- fmralign/piecewise.py | 467 +++++++++++++++++++++++++++++++ fmralign/tests/test_piecewise.py | 212 ++++++++++++++ 2 files changed, 679 insertions(+) create mode 100644 fmralign/piecewise.py create mode 100644 fmralign/tests/test_piecewise.py diff --git a/fmralign/piecewise.py b/fmralign/piecewise.py new file mode 100644 index 00000000..619ceb1a --- /dev/null +++ b/fmralign/piecewise.py @@ -0,0 +1,467 @@ +# -*- coding: utf-8 -*- +""" Test module to adapt quickly fmralign for piecewise srm +Implementation from fastSRM is taken from H. Richard +""" +# Author: T. Bazeille +# License: simplified BSD + +from nilearn.image import load_img +from sklearn.base import BaseEstimator, TransformerMixin +import numpy as np +from joblib import delayed, Parallel +from nilearn.image import concat_imgs, load_img +from nilearn.input_data.masker_validation import check_embedded_nifti_masker + +from ._utils import _make_parcellation, _intersect_clustering_mask +from sklearn.model_selection import ShuffleSplit +from sklearn.base import clone +import nibabel as nib +import warnings +import os + + +class Identity(BaseEstimator, TransformerMixin): + def fit(self, X, y=None): + self.basis_list = [np.eye(X[0].shape[0]) for _ in range(len(X))] + self.basis_list = [np.eye(X[0].shape[0]) for _ in range(len(X))] + return self + + def transform(self, X, y=None, subjects_indexes=None): + if subjects_indexes is None: + subjects_indexes = np.arange(len(self.basis_list)) + return np.array([X[i] for i in range(len(subjects_indexes))]) + + def inverse_transform(self, X, subjects_indexes=None): + return np.array([X for _ in subjects_indexes]) + + def add_subjects(self, X_list, S): + self.basis_list = [w for w in self.basis_list] + [ + np.eye(x.shape[0]) for x in X_list + ] + + +def generate_X_is(labels, X_list, masker, verbose): + """ Generate source and target data X_i and Y_i for each piece i. + + Parameters + ---------- + labels : list of ints (len n_features) + Parcellation of features in clusters + X_list: list of Niimgs + Source data + masker: instance of NiftiMasker or MultiNiftiMasker + Masker to be used on the data. For more information see: + http://nilearn.github.io/manipulating_images/masker_objects.html + verbose: integer, optional. + Indicate the level of verbosity. + + Yields + ------- + X_i: ndarray + Source data for piece i (shape : n_samples, n_features) + Y_i: ndarray + Target data for piece i (shape : n_samples, n_features) + + """ + + masked_X_list = [] + for X in X_list: + masked_X_list.append(masker.transform(X)) + unique_labels = np.unique(labels) + + for k in range(len(unique_labels)): + label = unique_labels[k] + i = label == labels + if (k + 1) % 25 == 0 and verbose > 0: + print( + "Fitting parcel: " + str(k + 1) + "/" + str(len(unique_labels)) + ) + # should return X_i Y_i + + yield [X_[:, i].T for X_ in masked_X_list] + + +def fit_one_piece(piece_X_list, method): + """ Align source and target data in one piece i, X_i and Y_i, using + alignment method and learn transformation to map X to Y. + + Parameters + ---------- + X_i: ndarray + Source data for piece i (shape : n_samples, n_features) + Y_i: ndarray + Target data for piece i (shape : n_samples, n_features) + alignment_method: string + Algorithm used to perform alignment between X_i and Y_i : + - either 'identity', 'scaled_orthogonal', 'optimal_transport', + 'ridge_cv', 'permutation', 'diagonal' + - or an instance of one of alignment classes + (imported from functional_alignment.alignment_methods) + Returns + ------- + alignment_algo + Instance of alignment estimator class fitted for X_i, Y_i + """ + + if method == "identity": + alignment_algo = Identity() + else: + from fastsrm.identifiable_srm import IdentifiableFastSRM + + # if isinstance(method, (FastSRM, MultiViewICA, AdaptiveMultiViewICA)): + if isinstance( + method, (IdentifiableFastSRM) + ): + alignment_algo = clone(method) + if hasattr(alignment_algo, "aggregate"): + alignment_algo.aggregate = None + if np.shape(piece_X_list)[1] < alignment_algo.n_components: + alignment_algo.n_components = np.shape(piece_X_list)[1] + else: + warn_msg = "Method not recognized, should be 'identity' or an instance of \ + FastSRM, MultiViewICA or AdaptiveMultiViewICA" + NotImplementedError(warn_msg) + # dirty monkey patching to avoid having n_components > n_voxels in any + #  piece which would yield a bug in add_subjects() + + reduced_SR = alignment_algo.fit(piece_X_list).transform(piece_X_list) + + if len(reduced_SR) == len(piece_X_list): + reduced_SR = np.mean(reduced_SR, axis=0) + + return alignment_algo, reduced_SR + + +def fit_one_parcellation( + X_list, + srm, + masker, + n_pieces, + clustering, + clustering_index, + n_jobs, + verbose, +): + """ Create one parcellation of n_pieces and align each source and target + data in one piece i, X_i and Y_i, using alignment method + and learn transformation to map X to Y. + + Parameters + ---------- + X_: Niimg-like object + Source data + alignment_method: string + algorithm used to perform alignment between each region of X_ and Y_ + masker: instance of NiftiMasker or MultiNiftiMasker + Masker to be used on the data. For more information see: + http://nilearn.github.io/manipulating_images/masker_objects.html + n_pieces: n_pieces: int, + Number of regions in which the data is parcellated for alignment + clustering: string or 3D Niimg + method used to perform parcellation of data. + If 3D Niimg, image used as predefined clustering. + clustering_index: list of integers + Clustering is performed on a 20% subset of the data chosen randomly + in timeframes. This index carry this subset. + n_jobs: integer, optional + The number of CPUs to use to do the computation. -1 means + 'all CPUs', -2 'all CPUs but one', and so on. + verbose: integer, optional + Indicate the level of verbosity. By default, nothing is printed + + Returns + ------- + alignment_algo + Instance of alignment estimator class fitted for X_i, Y_i + """ + # choose indexes maybe with index_img to not + labels = _make_parcellation( + X_list[0], + clustering_index, + clustering, + n_pieces, + masker, + verbose=verbose, + ) + + outputs = Parallel(n_jobs, prefer="threads", verbose=verbose)( + delayed(fit_one_piece)(piece_X_list, srm) + for piece_X_list in generate_X_is(labels, X_list, masker, verbose) + ) + fit = [output[0] for output in outputs] + reduced_sr = [output[1] for output in outputs] + return labels, fit, reduced_sr + + +class PiecewiseModel(BaseEstimator, TransformerMixin): + """ + Decompose the source images into regions and summarize subjects information \ + in a SR, then use alignment to predict \ + new contrast for target(s) subject. + """ + + def __init__( + self, + srm, + n_pieces=1, + clustering="kmeans", + n_bags=1, + mask=None, + smoothing_fwhm=None, + standardize=False, + detrend=None, + target_affine=None, + target_shape=None, + low_pass=None, + high_pass=None, + t_r=None, + n_jobs=1, + verbose=0, + reshape=True, + ): + """ + Parameters + ---------- + srm : FastSRM instance + n_pieces: int, optional (default = 1) + Number of regions in which the data is parcellated for alignment. + If 1 the alignment is done on full scale data. + If > 1, the voxels are clustered and alignment is performed on each + cluster applied to X and Y. + clustering : string or 3D Niimg optional (default : kmeans) + 'kmeans', 'ward', 'rena', 'hierarchical_kmeans' method used for + clustering of voxels based on functional signal, + passed to nilearn.regions.parcellations + If 3D Niimg, image used as predefined clustering, + n_bags and n_pieces are then ignored. + n_bags: int, optional (default = 1) + If 1 : one estimator is fitted. + If >1 number of bagged parcellations and estimators used. + mask: Niimg-like object, instance of NiftiMasker or + MultiNiftiMasker, optional (default = None) + Mask to be used on data. If an instance of masker is passed, then its + mask will be used. If no mask is given, it will be computed + automatically by a MultiNiftiMasker with default parameters. + smoothing_fwhm: float, optional (default = None) + If smoothing_fwhm is not None, it gives the size in millimeters + of the spatial smoothing to apply to the signal. + standardize: boolean, optional (default = None) + If standardize is True, the time-series are centered and normed: + their variance is put to 1 in the time dimension. + detrend: boolean, optional (default = None) + This parameter is passed to nilearn.signal.clean. + Please see the related documentation for details + target_affine: 3x3 or 4x4 matrix, optional (default = None) + This parameter is passed to nilearn.image.resample_img. + Please see the related documentation for details. + target_shape: 3-tuple of integers, optional (default = None) + This parameter is passed to nilearn.image.resample_img. + Please see the related documentation for details. + low_pass: None or float, optional (default = None) + This parameter is passed to nilearn.signal.clean. + Please see the related documentation for details. + high_pass: None or float, optional (default = None) + This parameter is passed to nilearn.signal.clean. + Please see the related documentation for details. + t_r: float, optional (default = None) + This parameter is passed to nilearn.signal.clean. + Please see the related documentation for details. + n_jobs: integer, optional (default = 1) + The number of CPUs to use to do the computation. -1 means + 'all CPUs', -2 'all CPUs but one', and so on. + verbose: integer, optional (default = 0) + Indicate the level of verbosity. By default, nothing is printed. + """ + self.srm = srm + self.n_pieces = n_pieces + self.clustering = clustering + self.n_bags = n_bags + self.mask = mask + self.smoothing_fwhm = smoothing_fwhm + self.standardize = standardize + self.detrend = detrend + self.target_affine = target_affine + self.target_shape = target_shape + self.low_pass = low_pass + self.high_pass = high_pass + self.t_r = t_r + self.n_jobs = n_jobs + self.memory = None + self.memory_level = 0 + self.verbose = verbose + self.reshape = reshape + + def fit(self, imgs): + """ + Learn a template from source images, using alignment. + + Parameters + ---------- + imgs: List of 4D Niimg-like or List of lists of 3D Niimg-like + Source subjects data. Each element of the parent list is one subject + data, and all must have the same length (n_samples). + + Returns + ------- + self + + Attributes + ---------- + self.template: 4D Niimg object + Length : n_samples + + """ + # Check if the input is a list, if list of lists, concatenate each subjects + # data into one unique image. + if not isinstance(imgs, (list, np.ndarray)) or len(imgs) < 2: + raise InputError( + "The method TemplateAlignment.fit() need a list input. \ + Each element of the list (Niimg-like or list of Niimgs) \ + is the data for one subject." + ) + else: + if isinstance(imgs[0], (list, np.ndarray)): + imgs = [concat_imgs(img) for img in imgs] + + self.masker_ = check_embedded_nifti_masker(self) + self.masker_.n_jobs = self.n_jobs # self.n_jobs + + # if masker_ has been provided a mask_img + if self.masker_.mask_img is None: + self.masker_.fit(imgs) + else: + self.masker_.fit() + + if type(self.clustering) == nib.nifti1.Nifti1Image or os.path.isfile( + self.clustering + ): + # check that clustering provided fills the mask, if not, reduce the mask + if 0 in self.masker_.transform(self.clustering): + reduced_mask = _intersect_clustering_mask( + self.clustering, self.masker_.mask_img + ) + self.mask = reduced_mask + self.masker_ = check_embedded_nifti_masker(self) + self.masker_.n_jobs = self.n_jobs + self.masker_.fit() + warnings.warn( + "Mask used was bigger than clustering provided. " + + "Its intersection with the clustering was used instead." + ) + + rs = ShuffleSplit(n_splits=self.n_bags, test_size=0.8, random_state=0) + + outputs = Parallel( + n_jobs=self.n_jobs, prefer="threads", verbose=self.verbose + )( + delayed(fit_one_parcellation)( + imgs, + self.srm, + self.masker_, + self.n_pieces, + self.clustering, + clustering_index, + self.n_jobs, + self.verbose, + ) + for clustering_index, _ in rs.split( + range(load_img(imgs[0]).shape[-1]) + ) + ) + + self.labels_ = [output[0] for output in outputs] + self.fit_ = [output[1] for output in outputs] + self.reduced_sr = [output[2] for output in outputs] + return self + + def add_subjects(self, imgs): + """ Add subject without recalculating SR""" + for labels, srm, reduced_sr in zip( + self.labels_, self.fit_, self.reduced_sr + ): + for X_i, piece_srm, piece_sr in zip( + list(generate_X_is(labels, imgs, self.masker_, self.verbose)), + srm, + reduced_sr, + ): + piece_srm.add_subjects(X_i, piece_sr) + return self + + def transform(self, imgs): + """ + imgs : list of Niimgs or string (paths). Masked shape : n_voxels, n_timeframes + + reshaped_aligned = + + !!!Not implemented for n_bags>1 + """ + n_comps = self.srm.n_components + aligned_imgs = [] + for labels, srm in zip(self.labels_, self.fit_): + bag_align = [] + for X_i, piece_srm in zip( + list(generate_X_is(labels, imgs, self.masker_, self.verbose)), + srm, + ): + piece_align = piece_srm.transform(X_i) + p_comps = piece_srm.n_components + if p_comps != n_comps: + piece_align = [ + np.pad( + t, + ((0, n_comps - p_comps), (0, 0)), + mode="constant", + ) + for t in piece_align + ] + bag_align.append(piece_align) + aligned_imgs.append(bag_align) + reordered_aligned = np.moveaxis(aligned_imgs[0], [1], [0]) + if self.reshape is False: + return reordered_aligned + reshaped_aligned = reordered_aligned.reshape( + len(imgs), -1, np.shape(aligned_imgs)[-1] + ) + return reshaped_aligned + + def get_full_basis(self): + """ + Concatenated (and padded if needed) local basis for each subject to + fullbrain basis of shape (n_components,n_voxels) for each subject + """ + unique_labels = np.unique(self.labels_[0]) + n_comps = self.srm.n_components + n_subs = len(self.fit_[0][0].basis_list) + full_basis_list = np.zeros( + shape=(n_subs, len(self.labels_[0]), n_comps) + ) + + for k in range(len(unique_labels)): + label = unique_labels[k] + k_estim = self.fit_[0][k] + i = label == self.labels_[0] + p_comps = k_estim.n_components + for s, basis in enumerate(k_estim.basis_list): + full_basis_list[s, i, :] = np.pad( + basis, ((0, n_comps - p_comps), (0, 0)), mode="constant" + ).T + + basis_imgs = [ + self.masker_.inverse_transform(full_basis.T) + for full_basis in full_basis_list + ] + return basis_imgs + + def fit_transform(self, imgs): + return self.fit(imgs).transform(imgs) + + def inverse_transform(self, X_transformed): + # TODO: Just inverse operations in transform to make API compliant + pass + + def clean(self): + """ + Just clean temporary directories + """ + for srm in self.fit_: + srm.clean() diff --git a/fmralign/tests/test_piecewise.py b/fmralign/tests/test_piecewise.py new file mode 100644 index 00000000..34c76b1e --- /dev/null +++ b/fmralign/tests/test_piecewise.py @@ -0,0 +1,212 @@ +from multiviewica import multiviewica +import pytest +import numpy as np +from sklearn.base import clone +from nilearn.input_data import NiftiMasker +from fmralign.piecewise import PiecewiseModel, Identity +from fastsrm.identifiable_srm import IdentifiableFastSRM + + +def to_niimgs(X, dim): + from nilearn.masking import _unmask_from_to_3d_array + import nibabel + + p = np.prod(dim) + assert len(dim) == 3 + assert X.shape[-1] <= p + mask = np.zeros(p).astype(np.bool) + mask[: X.shape[-1]] = 1 + assert mask.sum() == X.shape[1] + mask = mask.reshape(dim) + X = np.rollaxis( + np.array([_unmask_from_to_3d_array(x, mask) for x in X]), 0, start=4 + ) + affine = np.eye(4) + return ( + nibabel.Nifti1Image(X, affine), + nibabel.Nifti1Image(mask.astype(np.float), affine), + ) + + +n_bags = 1 +n_pieces = 1 +n_timeframes_align = 1000 +n_timeframes_test = 200 +n_subjects = 3 +# n_voxels should have an integer cubic root (8 ; 27 ; 64 ; 125 ...) +n_voxels = 8 + +grid_l = int(pow(n_voxels, 1.0 / 3.0)) + +_, mask = to_niimgs( + np.random.rand(n_timeframes_align, n_voxels), (grid_l, grid_l, grid_l) +) +masker = NiftiMasker(mask_img=mask).fit() + +train_data, test_data = {}, {} +for i in range(n_subjects): + train_data[i] = masker.inverse_transform( + np.random.rand(n_timeframes_align, 8) + ) + test_data[i] = masker.inverse_transform( + np.random.rand(n_timeframes_test, 8) + ) +masker = NiftiMasker(mask_img=mask).fit() + +n_components = 3 + + +srm = IdentifiableFastSRM(n_components=n_components, aggregate=None) +algos_to_test = [ + # "identity", + IdentifiableFastSRM(n_components=n_components, aggregate=None), +] + +@pytest.mark.parametrize("algo", algos_to_test) +def test_output_no_clustering(algo): + if algo == "identity": + n_comp = n_voxels + psrm = PiecewiseModel( + "identity", + n_pieces=n_pieces, + n_bags=n_bags, + clustering="kmeans", + mask=masker, + n_jobs=-1, + ) + algo = Identity() + else: + n_comp = algo.n_components + psrm = PiecewiseModel( + algo, + n_pieces=n_pieces, + n_bags=n_bags, + clustering="kmeans", + mask=masker, + n_jobs=-1, + ) + psrm.fit(list(train_data.values())[:-1]) + + srm_SR = algo.fit_transform( + [masker.transform(x).T for x in list(train_data.values())[:-1]] + ) + if len(srm_SR) == (n_subjects - 1): + srm_SR = np.mean(srm_SR, axis=0) + + np.shape(psrm.reduced_sr) + assert np.shape(psrm.reduced_sr) == ( + n_bags, + n_pieces, + n_comp, + n_timeframes_align, + ) + assert np.shape(psrm.labels_) == (n_bags, n_voxels) + assert np.shape(psrm.fit_) == (n_bags, n_pieces) + np.testing.assert_almost_equal(psrm.reduced_sr[0][0], srm_SR) + + algo.add_subjects( + [masker.transform(list(train_data.values())[-1]).T], srm_SR + ) + psrm.add_subjects([list(train_data.values())[-1]]) + + np.testing.assert_almost_equal(psrm.reduced_sr[0][0], srm_SR) + np.testing.assert_almost_equal(psrm.fit_[0][0].basis_list, algo.basis_list) + + aligned_test = psrm.transform(test_data.values()) + if hasattr(algo, "aggregate"): + algo.aggregate = None + srm_aligned_test = algo.transform( + [masker.transform(y).T for y in list(test_data.values())] + ) + assert np.shape(aligned_test) == (n_subjects, n_comp, n_timeframes_test) + np.testing.assert_almost_equal(aligned_test, srm_aligned_test) + + +def test_identity(): + id = Identity() + id_SR = id.fit_transform( + [masker.transform(x).T for x in list(train_data.values())[:-1]] + ) + id.add_subjects([masker.transform(list(train_data.values())[-1]).T], id_SR) + srm_aligned_test = id.transform( + [masker.transform(y).T for y in list(test_data.values())] + ) + # Check identity SRM just returns the data + np.testing.assert_almost_equal( + id_SR.shape, + np.asarray( + [masker.transform(x).T for x in list(train_data.values())[:-1]] + ).shape, + ) + + +@pytest.mark.parametrize("algo", algos_to_test) +def test_algo_each_piece(algo): + # Test that doing stuff piecewise give the same + # result as doing it separately + X = np.random.rand(1000, 8) + X1, mask = to_niimgs(X, (2, 2, 2)) + masker = NiftiMasker(mask_img=mask).fit() + X2 = masker.inverse_transform(np.random.rand(1000, 8)) + X3 = masker.inverse_transform(np.random.rand(1000, 8)) + X4 = masker.inverse_transform(np.random.rand(1000, 8)) + + cluster = np.array([1, 1, 1, 1, 2, 2, 2, 2]) + niimg_cluster = masker.inverse_transform(cluster) + if algo == "identity": + srm = PiecewiseModel("identity", mask=masker, clustering=niimg_cluster) + algo = Identity() + else: + srm = PiecewiseModel(algo, clustering=niimg_cluster, mask=masker) + S1 = np.array( + [ + clone(algo).fit_transform( + [masker.transform(x).T[cluster == i] for x in [X1, X2, X3, X4]] + ) + for i in [1, 2] + ] + ) + + S2 = srm.fit([X1, X2, X3, X4]).transform([X1, X2, X3, X4]) + S1 = np.swapaxes(S1, 0, 1) + S1 = S1.reshape(4, -1, 1000) + np.testing.assert_almost_equal(S1, S2) + + +n_components = 4 + + +@pytest.mark.parametrize("algo", algos_to_test) +def test_wrongshape(algo): + # Test that doing stuff piecewise give the same + # result as doing it separately + X = np.random.rand(10, 64) + X1, mask = to_niimgs(X, (4, 4, 4)) + masker = NiftiMasker(mask_img=mask).fit() + X2 = masker.inverse_transform(np.random.rand(10, 64)) + X3 = masker.inverse_transform(np.random.rand(10, 64)) + X4 = masker.inverse_transform(np.random.rand(10, 64)) + + cluster = np.array([1] * 8 + [2] * 56) + niimg_cluster = masker.inverse_transform(cluster) + if algo == "identity": + srm = PiecewiseModel("identity", mask=masker, clustering=niimg_cluster) + algo = Identity() + else: + srm = PiecewiseModel(algo, clustering=niimg_cluster, mask=masker) + S1 = np.array( + [ + clone(algo).fit_transform( + [ + masker.transform(x).T[cluster == i] + for x in [X1, X2, X3, X4] + ] + ) + for i in [1, 2] + ] + ) + + S2 = srm.fit([X1, X2, X3, X4]).transform([X1, X2, X3, X4]) + S1 = np.swapaxes(S1, 0, 1) + S1 = S1.reshape(4, -1, 10) + np.testing.assert_almost_equal(S1, S2) From b8990cc22204591399b731460375b99254b38527 Mon Sep 17 00:00:00 2001 From: hrichard Date: Mon, 29 Nov 2021 15:31:38 +0100 Subject: [PATCH 2/2] update dependencies --- README.md | 1 + fmralign/version.py | 4 ++++ 2 files changed, 5 insertions(+) diff --git a/README.md b/README.md index e4d1a3b2..4155467b 100644 --- a/README.md +++ b/README.md @@ -19,6 +19,7 @@ fmralign requires a Python installation and the following dependencies: * Matplotlib >= 3.1.1 * Nibabel >= 2.5.0 * Nilearn >= 1.5.0 +* FastSRM >= 0.0.4 ### Installation diff --git a/fmralign/version.py b/fmralign/version.py index f161a92a..9151bc0f 100644 --- a/fmralign/version.py +++ b/fmralign/version.py @@ -62,6 +62,10 @@ 'min_version': '0.5.0', 'required_at_installation': False, 'install_info': _FMRALIGN_INSTALL_MSG}) + ('fastsrm', { + 'min_version': '0.0.4', + 'required_at_installation': True, + 'install_info': _FMRALIGN_INSTALL_MSG}) )