diff --git a/pyproject.toml b/pyproject.toml index 6faf5185..8d876958 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,6 +28,7 @@ dependencies = [ "numpy>=1.17.3", "nest-asyncio>=1.5.1", "scikit-image>=0.14.2", + "scikit_learn", "scipy>=1.8.0", ] dynamic = ["version"] diff --git a/src/eddymotion/estimator.py b/src/eddymotion/estimator.py index 5c28f113..d274cfd8 100644 --- a/src/eddymotion/estimator.py +++ b/src/eddymotion/estimator.py @@ -120,6 +120,7 @@ def estimate( "avg", "average", "mean", + "gp", ) or model.lower().startswith("full") dwmodel = None diff --git a/src/eddymotion/model/__init__.py b/src/eddymotion/model/__init__.py index 3c44e2ad..ae454d40 100644 --- a/src/eddymotion/model/__init__.py +++ b/src/eddymotion/model/__init__.py @@ -26,6 +26,7 @@ AverageDWModel, DKIModel, DTIModel, + GPModel, ModelFactory, PETModel, TrivialB0Model, @@ -36,6 +37,7 @@ "AverageDWModel", "DKIModel", "DTIModel", + "GPModel", "TrivialB0Model", "PETModel", ) diff --git a/src/eddymotion/model/base.py b/src/eddymotion/model/base.py index 828d4496..4da06061 100644 --- a/src/eddymotion/model/base.py +++ b/src/eddymotion/model/base.py @@ -515,6 +515,13 @@ class DKIModel(BaseDWIModel): _model_class = "dipy.reconst.dki.DiffusionKurtosisModel" +class GPModel(BaseDWIModel): + """A wrapper of :obj:`~eddymotion.model.dipy.GaussianProcessModel`.""" + + _modelargs = ("kernel_model",) + _model_class = "eddymotion.model.dipy.GaussianProcessModel" + + def _rasb2dipy(gradient): gradient = np.asanyarray(gradient) if gradient.ndim == 1: diff --git a/src/eddymotion/model/dipy.py b/src/eddymotion/model/dipy.py new file mode 100644 index 00000000..1907363c --- /dev/null +++ b/src/eddymotion/model/dipy.py @@ -0,0 +1,419 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +# +# Copyright 2024 The NiPreps Developers +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# We support and encourage derived works from this project, please read +# about our expectations at +# +# https://www.nipreps.org/community/licensing/ +# +"""DIPY-like models (a sandbox to trial them out before upstreaming to DIPY).""" + +from __future__ import annotations + +import numpy as np +from dipy.core.gradients import GradientTable +from dipy.reconst.base import ReconstModel +from sklearn.gaussian_process import GaussianProcessRegressor + +# from dipy.reconst.multi_voxel import multi_voxel_fit + + +def gp_prediction( + model_gtab: GradientTable, + gtab: GradientTable, + model: GaussianProcessRegressor, + mask: np.ndarray | None = None, +) -> np.ndarray: + """ + Predicts one or more DWI orientations given a model. + + This function checks if the model is fitted and then extracts + orientations and potentially b-values from the gtab. It predicts the mean + and standard deviation of the DWI signal using the model. + + Parameters + ---------- + model: :obj:`~sklearn.gaussian_process.GaussianProcessRegressor` + A fitted GaussianProcessRegressor model. + gtab: :obj:`~dipy.core.gradients.GradientTable` + Gradient table with one or more orientations at which the GP will be evaluated. + mask: :obj:`numpy.ndarray` + A boolean mask indicating which voxels to use (optional). + + Returns + ------- + :obj:`numpy.ndarray` + A 3D or 4D array with the simulated gradient(s). + + """ + + # Check it's fitted as they do in sklearn internally + # https://github.com/scikit-learn/scikit-learn/blob/972e17fe1aa12d481b120ad4a3dc076bae736931/\ + # sklearn/gaussian_process/_gpr.py#L410C9-L410C42 + if not hasattr(model, "X_train_"): + raise RuntimeError("Model is not yet fitted.") + + # Extract orientations from gtab, and highly likely, the b-value too. + return model.predict(gtab, return_std=False) + + +def get_kernel( + kernel_model: str, + gtab: GradientTable | None = None, +) -> GaussianProcessRegressor.kernel: + """ + Returns a Gaussian process kernel based on the provided string. + + Currently supports 'test' kernel which is a combination of DotProduct and WhiteKernel + from scikit-learn. Raises a TypeError for unknown kernel models. + + Parameters + ---------- + kernel_model: :obj:`str` + The string representing the desired kernel model. + + Returns + ------- + :obj:`GaussianProcessRegressor.kernel` + A GaussianProcessRegressor kernel object. + + Raises + ------ + TypeError: If the provided ``kernel_model`` is not supported. + + """ + + if kernel_model == "spherical": + raise NotImplementedError("Spherical kernel is not currently implemented.") + + if kernel_model == "exponential": + raise NotImplementedError("Exponential kernel is not currently implemented.") + + if kernel_model == "test": + from sklearn.gaussian_process.kernels import DotProduct, WhiteKernel + + return DotProduct() + WhiteKernel() + + raise TypeError(f"Unknown kernel '{kernel_model}'.") + + +class GaussianProcessModel(ReconstModel): + """A Gaussian Process (GP) model to simulate single- and multi-shell DWI data.""" + + __slots__ = ( + "kernel_model", + "_modelfit", + ) + + def __init__( + self, + kernel_model: str = "spherical", + *args, + **kwargs, + ) -> None: + """A GP-based DWI model [Andersson15]_. + + Parameters + ---------- + kernel_model : :obj:`str` + Kernel model to calculate the GP's covariance matrix. + + References + ---------- + .. [Andersson15] Jesper L.R. Andersson and Stamatios N. Sotiropoulos. + Non-parametric representation and prediction of single- and multi-shell + diffusion-weighted MRI data using Gaussian processes. NeuroImage, 122:166-176, 2015. + doi:\ + `10.1016/j.neuroimage.2015.07.067 `__. + + """ + + ReconstModel.__init__(self, None) + self.kernel_model = kernel_model + + def fit( + self, + data: np.ndarray, + gtab: GradientTable | None = None, + mask: np.ndarray[bool] | None = None, + random_state: int = 0, + ) -> GPFit: + """Fit method of the DTI model class + + Parameters + ---------- + data : :obj:`~numpy.ndarray` + The measured signal from one voxel. + gtab : :obj:`~dipy.core.gradients.GradientTable` + The gradient table corresponding to the training data. + mask : :obj:`~numpy.ndarray` + A boolean array used to mark the coordinates in the data that + should be analyzed that has the shape data.shape[:-1] + + Returns + ------- + :obj:`~eddymotion.model.dipy.GPFit` + A model fit container object. + + """ + + data = ( + data[mask[..., None]] if mask is not None else np.reshape(data, (-1, data.shape[-1])) + ) + + signal_dirs = data.shape[-1] + grad_dirs = gtab.gradients.shape[0] + if signal_dirs != grad_dirs: + raise ValueError( + f"Mismatched data {signal_dirs} and gradient table {grad_dirs} sizes." + ) + + gpr = GaussianProcessRegressor( + kernel=get_kernel(self.kernel_model, gtab=gtab), + random_state=random_state, + ) + self._modelfit = GPFit( + gtab=gtab, + model=gpr.fit(gtab.gradients, data[0]), + mask=mask, + ) + return self._modelfit + + # @multi_voxel_fit + # def multi_fit(self, data_thres, mask=None, **kwargs): + # return GPFit(gpr.fit(self.gtab, data_thres)) + + def predict( + self, + gtab: GradientTable, + **kwargs, + ) -> np.ndarray: + """ + Predict using the Gaussian process model of the DWI signal for one or more gradients. + + Parameters + ---------- + gtab : :obj:`~dipy.core.gradients.GradientTable` + One or more gradient orientations at which the GP will be evaluated. + + Returns + ------- + :obj:`numpy.ndarray` + A 3D or 4D array with the simulated gradient(s). + + """ + return self._modelfit.predict(gtab) + + +class GPFit: + """ + A container class to store the fitted Gaussian process model and mask information. + + This class is typically returned by the `fit` and `multi_fit` methods of the + `GaussianProcessModel` class. It holds the fitted model and the mask used during fitting. + + Attributes + ---------- + fitted_gtab : :obj:`~dipy.core.gradients.GradientTable` + The original gradient table with which the GP has been fitted. + model: :obj:`~sklearn.gaussian_process.GaussianProcessRegressor` + The fitted Gaussian process regressor object. + mask: :obj:`~numpy.ndarray` + The boolean mask used during fitting (can be ``None``). + + """ + + def __init__( + self, + gtab: GradientTable, + model: GaussianProcessRegressor, + mask: np.ndarray | None = None, + ) -> None: + """ + Initialize a Gaussian Process fit container. + + Parameters + ---------- + gtab : :obj:`~dipy.core.gradients.GradientTable` + The gradient table with which the GP has been fitted. + model: :obj:`~sklearn.gaussian_process.GaussianProcessRegressor` + The fitted Gaussian process regressor object. + mask: :obj:`~numpy.ndarray` + The boolean mask used during fitting (can be ``None``). + + """ + self.fitted_gtab = gtab + self.model = model + self.mask = mask + + def predict( + self, + gtab: GradientTable, + ) -> np.ndarray: + """ + Generate DWI signal based on a fitted Gaussian Process. + + Parameters + ---------- + gtab: :obj:`~dipy.core.gradients.GradientTable` + Gradient table with one or more orientations at which the GP will be evaluated. + + Returns + ------- + :obj:`numpy.ndarray` + A 3D or 4D array with the simulated gradient(s). + + """ + return gp_prediction(self.fitted_gtab, gtab, self.model, mask=self.mask) + + +def _ensure_positive_scale( + a: float, +) -> None: + if a <= 0: + raise ValueError(f"a must be strictly positive. Provided: {a}") + + +def compute_exponential_covariance( + theta: np.ndarray, + a: float, +) -> np.ndarray: + r"""Compute the exponential covariance matrix following eq. 9 in [Andersson15]_. + + .. math:: + + C(\theta) = \exp(- \frac{\theta}{a}) + + Parameters + ---------- + theta : :obj:`~numpy.ndarray` + Pairwise angles across diffusion gradient encoding directions. + a : :obj:`float` + Positive scale parameter that here determines the "distance" at which θ + the covariance goes to zero. + + Returns + ------- + :obj:`~numpy.ndarray` + Exponential covariance function values. + + """ + + _ensure_positive_scale(a) + + return np.exp(-theta / a) + + +def compute_spherical_covariance( + theta: np.ndarray, + a: float, +) -> np.ndarray: + r"""Compute the spherical covariance matrix following eq. 10 in [Andersson15]_. + + .. math:: + + C(\theta) = \begin{cases} + 1 - \frac{3 \theta}{2 a} + \frac{\theta^3}{2 a^3} & \textnormal{if} \; \theta \leq a \\ + 0 & \textnormal{if} \; \theta > a + \end{cases} + + Parameters + ---------- + theta : :obj:`~numpy.ndarray` + Pairwise angles across diffusion gradient encoding directions. + a : :obj:`float` + Positive scale parameter that here determines the "distance" at which θ + the covariance goes to zero. + + Returns + ------- + :obj:`~numpy.ndarray` + Spherical covariance matrix. + + """ + + _ensure_positive_scale(a) + + return np.where(theta <= a, 1 - 3 * (theta / a) ** 2 + 2 * (theta / a) ** 3, 0) + + +def compute_pairwise_angles( + gtab_X: GradientTable, + gtab_Y: GradientTable | None = None, + closest_polarity: bool = True, +) -> np.ndarray: + r"""Compute pairwise angles across diffusion gradient encoding directions. + + Following [Andersson15]_, it computes the smallest of the angles between + each pair if ``closest_polarity`` is ``True``, i.e. + + .. math:: + + \theta(\mathbf{g}, \mathbf{g'}) = \arccos(\abs{\langle \mathbf{g}, \mathbf{g'} \rangle}) + + Parameters + ---------- + gtab_X: :obj:`~dipy.core.gradients.GradientTable` + Gradient table + gtab_Y: :obj:`~dipy.core.gradients.GradientTable` + Gradient table + closest_polarity : :obj:`bool` + ``True`` to consider the smallest of the two angles between the crossing + lines resulting from reversing each vector pair. + + Returns + ------- + :obj:`~numpy.ndarray` + Pairwise angles across diffusion gradient encoding directions. + + Examples + -------- + >>> from dipy.core.gradients import gradient_table + >>> bvecs = np.asarray([(1.0, -1.0), (0.0, 0.0), (0.0, 0.0)]) + >>> gtab = gradient_table([1000] * bvecs.shape[-1], bvecs) + >>> compute_pairwise_angles(gtab, closest_polarity=False)[0, 1] # doctest: +ELLIPSIS + 3.1415... + >>> bvecs1 = np.asarray([(1.0, -1.0), (0.0, 0.0), (0.0, 0.0)]) + >>> bvecs2 = np.asarray([(1.0, -1.0), (0.0, 0.0), (0.0, 0.0)]) + >>> gtab1 = gradient_table([1000] * bvecs1.shape[-1], bvecs1) + >>> gtab2 = gradient_table([1000] * bvecs2.shape[-1], bvecs2) + >>> compute_pairwise_angles(gtab1, gtab2, closest_polarity=False)[0, 1] # doctest: +ELLIPSIS + 3.1415... + >>> bvecs = np.asarray([(1.0, -1.0), (0.0, 0.0), (0.0, 0.0)]) + >>> gtab = gradient_table([1000] * bvecs.shape[-1], bvecs) + >>> compute_pairwise_angles(gtab)[0, 1] + 0.0 + + References + ---------- + .. [Andersson15] J. L. R. Andersson. et al., An integrated approach to + correction for off-resonance effects and subject movement in diffusion MR + imaging, NeuroImage 125 (2016) 1063–1078 + + """ + + bvecs_X = gtab_X.bvecs.T + bvecs_X = np.array(bvecs_X) / np.linalg.norm(bvecs_X, axis=0) + + if gtab_Y is None: + bvecs_Y = bvecs_X + else: + bvecs_Y = gtab_Y.bvecs.T + bvecs_Y = np.array(bvecs_Y) / np.linalg.norm(bvecs_Y, axis=0) + + cosines = np.clip(bvecs_X.T @ bvecs_Y, -1.0, 1.0) + return np.arccos(np.abs(cosines) if closest_polarity else cosines) diff --git a/src/eddymotion/model/dmri_covariance.py b/src/eddymotion/model/dmri_covariance.py deleted file mode 100644 index 4af32e84..00000000 --- a/src/eddymotion/model/dmri_covariance.py +++ /dev/null @@ -1,83 +0,0 @@ -# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- -# vi: set ft=python sts=4 ts=4 sw=4 et: -# -# Copyright 2024 The NiPreps Developers -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# We support and encourage derived works from this project, please read -# about our expectations at -# -# https://www.nipreps.org/community/licensing/ -# -import numpy as np - - -def _ensure_positive_scale(a): - if a <= 0: - raise ValueError(f"a must be strictly positive. Provided: {a}") - - -def compute_exponential_covariance(theta, a): - r"""Compute the exponential covariance matrix following eq. 9 in [Andersson15]_. - - .. math:: - - C(\theta) = \exp(- \frac{\theta}{a}) - - Parameters - ---------- - theta : :obj:`~numpy.ndarray` - Pairwise angles across diffusion gradient encoding directions. - a : :obj:`float` - Positive scale parameter that here determines the "distance" at which θ - the covariance goes to zero. - - Returns - ------- - :obj:`~numpy.ndarray` - Exponential covariance function values. - """ - - _ensure_positive_scale(a) - - return np.exp(-theta / a) - - -def compute_spherical_covariance(theta, a): - r"""Compute the spherical covariance matrix following eq. 10 in [Andersson15]_. - - .. math:: - - C(\theta) = \begin{cases} - 1 - \frac{3 \theta}{2 a} + \frac{\theta^3}{2 a^3} & \textnormal{if} \; \theta \leq a \\ - 0 & \textnormal{if} \; \theta > a - \end{cases} - - Parameters - ---------- - theta : :obj:`~numpy.ndarray` - Pairwise angles across diffusion gradient encoding directions. - a : :obj:`float` - Positive scale parameter that here determines the "distance" at which θ - the covariance goes to zero. - - Returns - ------- - :obj:`~numpy.ndarray` - Spherical covariance matrix. - """ - - _ensure_positive_scale(a) - - return np.where(theta <= a, 1 - 3 * (theta / a) ** 2 + 2 * (theta / a) ** 3, 0) diff --git a/src/eddymotion/model/gradient_utils.py b/src/eddymotion/model/gradient_utils.py deleted file mode 100644 index 1d93daa3..00000000 --- a/src/eddymotion/model/gradient_utils.py +++ /dev/null @@ -1,83 +0,0 @@ -# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- -# vi: set ft=python sts=4 ts=4 sw=4 et: -# -# Copyright 2024 The NiPreps Developers -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# We support and encourage derived works from this project, please read -# about our expectations at -# -# https://www.nipreps.org/community/licensing/ -# -import numpy as np - - -def compute_pairwise_angles(bvecs1, bvecs2, closest_polarity): - r"""Compute pairwise angles across diffusion gradient encoding directions. - - Following [Andersson15]_, it computes the smallest of the angles between - each pair if ``closest_polarity`` is ``True``, i.e. - - .. math:: - - \theta(\mathbf{g}, \mathbf{g'}) = \arccos(\abs{\langle \mathbf{g}, \mathbf{g'} \rangle}) - - Parameters - ---------- - bvecs1 : :obj:`~numpy.ndarray` - Diffusion gradient encoding directions in FSL format. - bvecs2 : :obj:`~numpy.ndarray` - Diffusion gradient encoding directions in FSL format. - closest_polarity : :obj:`bool` - ``True`` to consider the smallest of the two angles between the crossing - lines resulting from reversing each vector pair. - - Returns - ------- - :obj:`~numpy.ndarray` - Pairwise angles across diffusion gradient encoding directions. - - Examples - -------- - >>> compute_pairwise_angles( - ... ((1.0, -1.0), (0.0, 0.0), (0.0, 0.0)), - ... ((1.0, -1.0), (0.0, 0.0), (0.0, 0.0)), - ... False, - ... )[0, 1] # doctest: +ELLIPSIS - 3.1415... - >>> compute_pairwise_angles( - ... ((1.0, -1.0), (0.0, 0.0), (0.0, 0.0)), - ... ((1.0, -1.0), (0.0, 0.0), (0.0, 0.0)), - ... True, - ... )[0, 1] - 0.0 - - References - ---------- - .. [Andersson15] J. L. R. Andersson. et al., An integrated approach to - correction for off-resonance effects and subject movement in diffusion MR - imaging, NeuroImage 125 (2016) 1063–1078 - """ - - if np.shape(bvecs1)[0] != 3: - raise ValueError(f"bvecs1 must be of shape (3, N). Found: {bvecs1.shape}") - - if np.shape(bvecs2)[0] != 3: - raise ValueError(f"bvecs2 must be of shape (3, N). Found: {bvecs2.shape}") - - # Ensure b-vectors are unit-norm - bvecs1 = np.array(bvecs1) / np.linalg.norm(bvecs1, axis=0) - bvecs2 = np.array(bvecs2) / np.linalg.norm(bvecs2, axis=0) - cosines = np.clip(bvecs1.T @ bvecs2, -1.0, 1.0) - return np.arccos(np.abs(cosines) if closest_polarity else cosines) diff --git a/src/eddymotion/model/tests/test_dmri_covariance.py b/src/eddymotion/model/tests/test_dmri_covariance.py deleted file mode 100644 index 0d75a1fc..00000000 --- a/src/eddymotion/model/tests/test_dmri_covariance.py +++ /dev/null @@ -1,65 +0,0 @@ -# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- -# vi: set ft=python sts=4 ts=4 sw=4 et: -# -# Copyright 2024 The NiPreps Developers -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# We support and encourage derived works from this project, please read -# about our expectations at -# -# https://www.nipreps.org/community/licensing/ -# -import numpy as np -import pytest - -from eddymotion.model.dmri_covariance import ( - compute_exponential_covariance, - compute_spherical_covariance, -) - - -@pytest.mark.parametrize( - ("theta", "a", "expected"), - [ - ( - np.asarray( - [0.0, np.pi / 2, np.pi / 2, np.pi / 4, np.pi / 4, np.pi / 2, np.pi / 4], - ), - 1.0, - np.asarray( - [1.0, 0.20787958, 0.20787958, 0.45593813, 0.45593813, 0.20787958, 0.45593813] - ), - ) - ], -) -def test_compute_exponential_covariance(theta, a, expected): - obtained = compute_exponential_covariance(theta, a) - assert np.allclose(obtained, expected) - - -@pytest.mark.parametrize( - ("theta", "a", "expected"), - [ - ( - np.asarray( - [0.0, np.pi / 2, np.pi / 2, np.pi / 4, np.pi / 4, np.pi / 2, np.pi / 4], - ), - 1.0, - np.asarray([1.0, 0.0, 0.0, 0.11839532, 0.11839532, 0.0, 0.11839532]), - ) - ], -) -def test_compute_spherical_covariance(theta, a, expected): - obtained = compute_spherical_covariance(theta, a) - assert np.allclose(obtained, expected) diff --git a/src/eddymotion/model/tests/test_gradient_utils.py b/test/test_dipy.py similarity index 60% rename from src/eddymotion/model/tests/test_gradient_utils.py rename to test/test_dipy.py index b7be16e7..07b257f9 100644 --- a/src/eddymotion/model/tests/test_gradient_utils.py +++ b/test/test_dipy.py @@ -22,8 +22,13 @@ # import numpy as np import pytest +from dipy.core.gradients import gradient_table -from eddymotion.model.gradient_utils import compute_pairwise_angles +from eddymotion.model.dipy import ( + compute_exponential_covariance, + compute_pairwise_angles, + compute_spherical_covariance, +) # No need to use normalized vectors: compute_pairwise_angles takes care of it. @@ -41,6 +46,21 @@ [0, 0, 1, 0, 1, 1, 1], ] ), + None, + True, + np.array( + [ + [0.0, np.pi / 2, np.pi / 2, np.pi / 4, np.pi / 4, np.pi / 2, np.pi / 4], + [np.pi / 2, 0.0, np.pi / 2, np.pi / 4, np.pi / 2, np.pi / 4, np.pi / 2], + [np.pi / 2, np.pi / 2, 0.0, np.pi / 2, np.pi / 4, np.pi / 4, np.pi / 4], + [np.pi / 4, np.pi / 4, np.pi / 2, 0.0, np.pi / 3, np.pi / 3, np.pi / 3], + [np.pi / 4, np.pi / 2, np.pi / 4, np.pi / 3, 0.0, np.pi / 3, np.pi / 2], + [np.pi / 2, np.pi / 4, np.pi / 4, np.pi / 3, np.pi / 3, 0.0, np.pi / 3], + [np.pi / 4, np.pi / 2, np.pi / 4, np.pi / 3, np.pi / 2, np.pi / 3, 0.0], + ] + ), + ), + ( np.array( [ [1, 0, 0, 1, 1, 0, -1], @@ -48,16 +68,25 @@ [0, 0, 1, 0, 1, 1, 1], ] ), - True, + None, + False, np.array( [ - [0.0, np.pi / 2, np.pi / 2, np.pi / 4, np.pi / 4, np.pi / 2, np.pi / 4], + [0.0, np.pi / 2, np.pi / 2, np.pi / 4, np.pi / 4, np.pi / 2, 3 * np.pi / 4], [np.pi / 2, 0.0, np.pi / 2, np.pi / 4, np.pi / 2, np.pi / 4, np.pi / 2], [np.pi / 2, np.pi / 2, 0.0, np.pi / 2, np.pi / 4, np.pi / 4, np.pi / 4], - [np.pi / 4, np.pi / 4, np.pi / 2, 0.0, np.pi / 3, np.pi / 3, np.pi / 3], + [np.pi / 4, np.pi / 4, np.pi / 2, 0.0, np.pi / 3, np.pi / 3, 2 * np.pi / 3], [np.pi / 4, np.pi / 2, np.pi / 4, np.pi / 3, 0.0, np.pi / 3, np.pi / 2], [np.pi / 2, np.pi / 4, np.pi / 4, np.pi / 3, np.pi / 3, 0.0, np.pi / 3], - [np.pi / 4, np.pi / 2, np.pi / 4, np.pi / 3, np.pi / 2, np.pi / 3, 0.0], + [ + 3 * np.pi / 4, + np.pi / 2, + np.pi / 4, + 2 * np.pi / 3, + np.pi / 2, + np.pi / 3, + 0.0, + ], ] ), ), @@ -120,8 +149,55 @@ ], ) def test_compute_pairwise_angles(bvecs1, bvecs2, closest_polarity, expected): - obtained = compute_pairwise_angles(bvecs1, bvecs2, closest_polarity) + # DIPY requires the vectors to be normalized + _bvecs1 = bvecs1 / np.linalg.norm(bvecs1, axis=0) + gtab1 = gradient_table([1000] * _bvecs1.shape[-1], _bvecs1) - assert (bvecs1.shape[-1], bvecs2.shape[-1]) == obtained.shape + _bvecs2 = None + gtab2 = None + if bvecs2 is not None: + _bvecs2 = bvecs2 / np.linalg.norm(bvecs2, axis=0) + gtab2 = gradient_table([1000] * _bvecs2.shape[-1], _bvecs2) + + obtained = compute_pairwise_angles(gtab1, gtab2, closest_polarity) + + if _bvecs2 is not None: + assert (_bvecs1.shape[-1], _bvecs2.shape[-1]) == obtained.shape assert obtained.shape == expected.shape np.testing.assert_array_almost_equal(obtained, expected, decimal=2) + + +@pytest.mark.parametrize( + ("theta", "a", "expected"), + [ + ( + np.asarray( + [0.0, np.pi / 2, np.pi / 2, np.pi / 4, np.pi / 4, np.pi / 2, np.pi / 4], + ), + 1.0, + np.asarray( + [1.0, 0.20787958, 0.20787958, 0.45593813, 0.45593813, 0.20787958, 0.45593813] + ), + ) + ], +) +def test_compute_exponential_covariance(theta, a, expected): + obtained = compute_exponential_covariance(theta, a) + assert np.allclose(obtained, expected) + + +@pytest.mark.parametrize( + ("theta", "a", "expected"), + [ + ( + np.asarray( + [0.0, np.pi / 2, np.pi / 2, np.pi / 4, np.pi / 4, np.pi / 2, np.pi / 4], + ), + 1.0, + np.asarray([1.0, 0.0, 0.0, 0.11839532, 0.11839532, 0.0, 0.11839532]), + ) + ], +) +def test_compute_spherical_covariance(theta, a, expected): + obtained = compute_spherical_covariance(theta, a) + assert np.allclose(obtained, expected) diff --git a/test/test_model.py b/test/test_model.py index 7c7a906f..dba71522 100644 --- a/test/test_model.py +++ b/test/test_model.py @@ -24,12 +24,15 @@ import numpy as np import pytest +from dipy.core.gradients import gradient_table +from sklearn.datasets import make_regression from eddymotion import model from eddymotion.data.dmri import DWI from eddymotion.data.splitting import lovo_split from eddymotion.exceptions import ModelNotFittedError from eddymotion.model.base import DEFAULT_MAX_S0, DEFAULT_MIN_S0 +from eddymotion.model.dipy import GaussianProcessModel def test_trivial_model(): @@ -105,6 +108,22 @@ def test_average_model(): assert np.all(tmodel_2000.predict([0, 0, 0]) == 1100) +def test_gp_model(): + gp = GaussianProcessModel("test") + + assert isinstance(gp, model.dipy.GaussianProcessModel) + + X, y = make_regression(n_samples=100, n_features=3, noise=0, random_state=0) + + bvecs = X.T / np.linalg.norm(X.T, axis=0) + gtab = gradient_table([1000] * bvecs.shape[-1], bvecs) + gp.fit(y, gtab) + X_qry = bvecs[:, :2].T + prediction = gp.predict(X_qry, return_std=True) + + assert prediction.shape == (X_qry.shape[0],) + + def test_two_initialisations(datadir): """Check that the two different initialisations result in the same models"""