From e90ebae30eb17ca92d7ea82be0a040a6a3c40c9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jon=20Haitz=20Legarreta=20Gorro=C3=B1o?= Date: Sun, 7 Apr 2024 10:36:04 -0400 Subject: [PATCH 1/9] ENH: Implement Gaussian Process Implement Gaussian Process. --- pyproject.toml | 1 + src/eddymotion/estimator.py | 1 + src/eddymotion/model/__init__.py | 2 + src/eddymotion/model/base.py | 67 ++++++++++++++++++++++++++++++++ test/test_model.py | 19 +++++++++ 5 files changed, 90 insertions(+) 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..b64712ff 100644 --- a/src/eddymotion/model/__init__.py +++ b/src/eddymotion/model/__init__.py @@ -26,6 +26,7 @@ AverageDWModel, DKIModel, DTIModel, + GaussianProcessModel, ModelFactory, PETModel, TrivialB0Model, @@ -36,6 +37,7 @@ "AverageDWModel", "DKIModel", "DTIModel", + "GaussianProcessModel", "TrivialB0Model", "PETModel", ) diff --git a/src/eddymotion/model/base.py b/src/eddymotion/model/base.py index 828d4496..69986112 100644 --- a/src/eddymotion/model/base.py +++ b/src/eddymotion/model/base.py @@ -27,6 +27,7 @@ import numpy as np from dipy.core.gradients import gradient_table from joblib import Parallel, delayed +from sklearn.gaussian_process import GaussianProcessRegressor from eddymotion.exceptions import ModelNotFittedError @@ -515,6 +516,72 @@ class DKIModel(BaseDWIModel): _model_class = "dipy.reconst.dki.DiffusionKurtosisModel" +class GaussianProcessModel(BaseModel): + """A Gaussian process model for DWI data based on [Andersson15]_.""" + + __slots__ = ( + "_dwi", + "_kernel", + "_gpr", + ) + + def __init__(self, dwi, kernel, **kwargs): + """Implement object initialization. + + Parameters + ---------- + dwi : :obj:`~eddymotion.dmri.DWI` + The DWI data. + kernel : :obj:`~sklearn.gaussian_process.kernels.Kernel` + Kernel instance. + """ + + self._dwi = dwi + self._kernel = kernel + + def fit(self, X, y, *args, **kwargs): + """Fit the Gaussian process model to the training data. + + Parameters + ---------- + X : :obj:`~numpy.ndarray` of shape (n_samples, n_features) + Feature values for training. For the DWI cae, ``n_samples`` is the + number of diffusion-encoding gradient vectors, and ``n_features`` + being 3 (the spatial coordinates). + y : :obj:`~numpy.ndarray` of shape (n_samples,) or (n_samples, n_targets) + Target values: the DWI signal values. + """ + + self._gpr = GaussianProcessRegressor(kernel=self._kernel, random_state=0) + self._gpr.fit(X, y) + self._is_fitted = True + + def predict(self, X, **kwargs): + """Predict using the Gaussian process model of the DWI signal, where + ``X`` is a diffusion-encoding gradient vector whose DWI data needs to be + estimated. + + Parameters + ---------- + X : :obj:`~numpy.ndarray` of shape (n_samples,) + Query points where the Gaussian process is evaluated: the + diffusion-encoding gradient vectors of interest. + + Returns + ------- + y_mean : :obj:`~numpy.ndarray` of shape (n_samples,) or (n_samples, n_targets) + Mean of predictive distribution at query points. + y_std : :obj:`~numpy.ndarray` of shape (n_samples,) or (n_samples, n_targets) + Standard deviation of predictive distribution at query points. + """ + + if not self._is_fitted: + raise ModelNotFittedError(f"{type(self).__name__} must be fitted before predicting") + + y_mean, y_std = self._gpr.predict(X, return_std=True) + return y_mean, y_std + + def _rasb2dipy(gradient): gradient = np.asanyarray(gradient) if gradient.ndim == 1: diff --git a/test/test_model.py b/test/test_model.py index 7c7a906f..23a097b0 100644 --- a/test/test_model.py +++ b/test/test_model.py @@ -24,6 +24,8 @@ import numpy as np import pytest +from sklearn.datasets import make_friedman2 +from sklearn.gaussian_process.kernels import DotProduct, WhiteKernel from eddymotion import model from eddymotion.data.dmri import DWI @@ -105,6 +107,23 @@ def test_average_model(): assert np.all(tmodel_2000.predict([0, 0, 0]) == 1100) +def test_gp_model(datadir): + dwi = DWI.from_filename(datadir / "dwi.h5") + + kernel = DotProduct() + WhiteKernel() + + gp = model.GaussianProcessModel(dwi=dwi, kernel=kernel) + + assert isinstance(gp, model.GaussianProcessModel) + + X, y = make_friedman2(n_samples=500, noise=0, random_state=0) + gp.fit(X, y) + X_qry = X[:2, :] + 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""" From 712eaef6b2b5824c1636e4976d73c342585d22b0 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 4 Jul 2024 15:49:37 +0200 Subject: [PATCH 2/9] ENH: Refactor the model so it is compatible with DIPY's interface (#7) Refactor the model so that it is compatible with DIPY's reconstruction interface. Bring utilities into a new `dipy` module so that every part required for the Gaussian process (gradient angle computation, covariance methods, etc.) is more localized without proliferating small modules. Co-authored-by: Ariel Rokem --- src/eddymotion/model/base.py | 68 +-- src/eddymotion/model/dipy.py | 417 ++++++++++++++++++ src/eddymotion/model/dmri_covariance.py | 83 ---- src/eddymotion/model/gradient_utils.py | 83 ---- .../model/tests/test_dmri_covariance.py | 65 --- .../test_dipy.py | 42 +- 6 files changed, 462 insertions(+), 296 deletions(-) create mode 100644 src/eddymotion/model/dipy.py delete mode 100644 src/eddymotion/model/dmri_covariance.py delete mode 100644 src/eddymotion/model/gradient_utils.py delete mode 100644 src/eddymotion/model/tests/test_dmri_covariance.py rename src/eddymotion/model/tests/test_gradient_utils.py => test/test_dipy.py (79%) diff --git a/src/eddymotion/model/base.py b/src/eddymotion/model/base.py index 69986112..4da06061 100644 --- a/src/eddymotion/model/base.py +++ b/src/eddymotion/model/base.py @@ -27,7 +27,6 @@ import numpy as np from dipy.core.gradients import gradient_table from joblib import Parallel, delayed -from sklearn.gaussian_process import GaussianProcessRegressor from eddymotion.exceptions import ModelNotFittedError @@ -516,70 +515,11 @@ class DKIModel(BaseDWIModel): _model_class = "dipy.reconst.dki.DiffusionKurtosisModel" -class GaussianProcessModel(BaseModel): - """A Gaussian process model for DWI data based on [Andersson15]_.""" +class GPModel(BaseDWIModel): + """A wrapper of :obj:`~eddymotion.model.dipy.GaussianProcessModel`.""" - __slots__ = ( - "_dwi", - "_kernel", - "_gpr", - ) - - def __init__(self, dwi, kernel, **kwargs): - """Implement object initialization. - - Parameters - ---------- - dwi : :obj:`~eddymotion.dmri.DWI` - The DWI data. - kernel : :obj:`~sklearn.gaussian_process.kernels.Kernel` - Kernel instance. - """ - - self._dwi = dwi - self._kernel = kernel - - def fit(self, X, y, *args, **kwargs): - """Fit the Gaussian process model to the training data. - - Parameters - ---------- - X : :obj:`~numpy.ndarray` of shape (n_samples, n_features) - Feature values for training. For the DWI cae, ``n_samples`` is the - number of diffusion-encoding gradient vectors, and ``n_features`` - being 3 (the spatial coordinates). - y : :obj:`~numpy.ndarray` of shape (n_samples,) or (n_samples, n_targets) - Target values: the DWI signal values. - """ - - self._gpr = GaussianProcessRegressor(kernel=self._kernel, random_state=0) - self._gpr.fit(X, y) - self._is_fitted = True - - def predict(self, X, **kwargs): - """Predict using the Gaussian process model of the DWI signal, where - ``X`` is a diffusion-encoding gradient vector whose DWI data needs to be - estimated. - - Parameters - ---------- - X : :obj:`~numpy.ndarray` of shape (n_samples,) - Query points where the Gaussian process is evaluated: the - diffusion-encoding gradient vectors of interest. - - Returns - ------- - y_mean : :obj:`~numpy.ndarray` of shape (n_samples,) or (n_samples, n_targets) - Mean of predictive distribution at query points. - y_std : :obj:`~numpy.ndarray` of shape (n_samples,) or (n_samples, n_targets) - Standard deviation of predictive distribution at query points. - """ - - if not self._is_fitted: - raise ModelNotFittedError(f"{type(self).__name__} must be fitted before predicting") - - y_mean, y_std = self._gpr.predict(X, return_std=True) - return y_mean, y_std + _modelargs = ("kernel_model",) + _model_class = "eddymotion.model.dipy.GaussianProcessModel" def _rasb2dipy(gradient): diff --git a/src/eddymotion/model/dipy.py b/src/eddymotion/model/dipy.py new file mode 100644 index 00000000..6495962c --- /dev/null +++ b/src/eddymotion/model/dipy.py @@ -0,0 +1,417 @@ +# 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._gpr, "X_train_"): + raise RuntimeError("Model is not yet fitted.") + + # Extract orientations from gtab, and highly likely, the b-value too. + return model._gpr.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])) + ) + + if data.shape[-1] != len(gtab): + raise ValueError( + f"Mismatched data {data.shape[-1]} and gradient table {len(gtab)} sizes." + ) + + gpr = GaussianProcessRegressor( + kernel=get_kernel(self.kernel_model, gtab=gtab), + random_state=random_state, + ) + self._modelfit = GPFit( + gpr.fit(gtab.gradients, data), + gtab=gtab, + 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 + -------- + >>> compute_pairwise_angles( + ... ((1.0, -1.0), (0.0, 0.0), (0.0, 0.0)), + ... closest_polarity=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)), + ... closest_polarity=False, + ... )[0, 1] # doctest: +ELLIPSIS + 3.1415... + >>> compute_pairwise_angles( + ... ((1.0, -1.0), (0.0, 0.0), (0.0, 0.0)), + ... )[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 79% rename from src/eddymotion/model/tests/test_gradient_utils.py rename to test/test_dipy.py index b7be16e7..a4bb79ff 100644 --- a/src/eddymotion/model/tests/test_gradient_utils.py +++ b/test/test_dipy.py @@ -23,7 +23,11 @@ import numpy as np import pytest -from eddymotion.model.gradient_utils import compute_pairwise_angles +from eddymotion.model.dipy import ( + compute_pairwise_angles, + compute_exponential_covariance, + compute_spherical_covariance, +) # No need to use normalized vectors: compute_pairwise_angles takes care of it. @@ -125,3 +129,39 @@ def test_compute_pairwise_angles(bvecs1, bvecs2, closest_polarity, expected): 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) From fb4f87c8234e05a9f4f2e54f47adc611ace17d27 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 4 Jul 2024 16:05:15 +0200 Subject: [PATCH 3/9] sty: fix import order --- test/test_dipy.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_dipy.py b/test/test_dipy.py index a4bb79ff..5ebebbf8 100644 --- a/test/test_dipy.py +++ b/test/test_dipy.py @@ -24,8 +24,8 @@ import pytest from eddymotion.model.dipy import ( - compute_pairwise_angles, compute_exponential_covariance, + compute_pairwise_angles, compute_spherical_covariance, ) From 832ca762efd32443269ed85748e7299aba5012b3 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 4 Jul 2024 16:07:11 +0200 Subject: [PATCH 4/9] fix: amend bad import --- src/eddymotion/model/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/eddymotion/model/__init__.py b/src/eddymotion/model/__init__.py index b64712ff..ae454d40 100644 --- a/src/eddymotion/model/__init__.py +++ b/src/eddymotion/model/__init__.py @@ -26,7 +26,7 @@ AverageDWModel, DKIModel, DTIModel, - GaussianProcessModel, + GPModel, ModelFactory, PETModel, TrivialB0Model, @@ -37,7 +37,7 @@ "AverageDWModel", "DKIModel", "DTIModel", - "GaussianProcessModel", + "GPModel", "TrivialB0Model", "PETModel", ) From a72b9abd5befaa116cf2c61a26d0e04360202739 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 4 Jul 2024 16:32:01 +0200 Subject: [PATCH 5/9] fix: meet new signature --- test/test_dipy.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/test/test_dipy.py b/test/test_dipy.py index 5ebebbf8..c25fedef 100644 --- a/test/test_dipy.py +++ b/test/test_dipy.py @@ -22,6 +22,7 @@ # import numpy as np import pytest +from dipy.core.gradients import gradient_table from eddymotion.model.dipy import ( compute_exponential_covariance, @@ -124,7 +125,9 @@ ], ) def test_compute_pairwise_angles(bvecs1, bvecs2, closest_polarity, expected): - obtained = compute_pairwise_angles(bvecs1, bvecs2, closest_polarity) + gtab1 = gradient_table([1000] * len(bvecs1), bvecs1) + gtab2 = gradient_table([1000] * len(bvecs2), bvecs2) + obtained = compute_pairwise_angles(gtab1, gtab2, closest_polarity) assert (bvecs1.shape[-1], bvecs2.shape[-1]) == obtained.shape assert obtained.shape == expected.shape From e75b6e56e746508706a3719b1ec48ffe6bcb371f Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 4 Jul 2024 16:32:20 +0200 Subject: [PATCH 6/9] fix: revise test of gaussian process --- test/test_model.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/test/test_model.py b/test/test_model.py index 23a097b0..9f1fdc6e 100644 --- a/test/test_model.py +++ b/test/test_model.py @@ -25,13 +25,13 @@ import numpy as np import pytest from sklearn.datasets import make_friedman2 -from sklearn.gaussian_process.kernels import DotProduct, WhiteKernel 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(): @@ -107,12 +107,8 @@ def test_average_model(): assert np.all(tmodel_2000.predict([0, 0, 0]) == 1100) -def test_gp_model(datadir): - dwi = DWI.from_filename(datadir / "dwi.h5") - - kernel = DotProduct() + WhiteKernel() - - gp = model.GaussianProcessModel(dwi=dwi, kernel=kernel) +def test_gp_model(): + gp = GaussianProcessModel(kernel="default") assert isinstance(gp, model.GaussianProcessModel) From 0c32def2f307f4f2bf815c0ab1f388e6e802b322 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jon=20Haitz=20Legarreta=20Gorro=C3=B1o?= Date: Thu, 4 Jul 2024 11:44:15 -0400 Subject: [PATCH 7/9] BUG: Fix tests Fix GP model testing: `GaussianProcessModel` is now hosted in the `model.dipy` module. Fix bvals array dimensionality and make bvecs be unit vectors in angle computation test. Fix the pairwise angle computation doctests. Test the case where the second argument to the pairwaise angle computation function is `None`. Simplify the test parameterization as the former corresponds to providing the same value for both gtab arguments. Adapt the test so that we avoid making computations if the second argument is `None`. --- src/eddymotion/model/dipy.py | 24 +++++++++--------- test/test_dipy.py | 47 ++++++++++++++++++++++++++++++------ test/test_model.py | 2 +- 3 files changed, 53 insertions(+), 20 deletions(-) diff --git a/src/eddymotion/model/dipy.py b/src/eddymotion/model/dipy.py index 6495962c..63e2150f 100644 --- a/src/eddymotion/model/dipy.py +++ b/src/eddymotion/model/dipy.py @@ -380,20 +380,20 @@ def compute_pairwise_angles( Examples -------- - >>> compute_pairwise_angles( - ... ((1.0, -1.0), (0.0, 0.0), (0.0, 0.0)), - ... closest_polarity=False, - ... )[0, 1] # doctest: +ELLIPSIS + >>> 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... - >>> 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)), - ... closest_polarity=False, - ... )[0, 1] # doctest: +ELLIPSIS + >>> 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... - >>> compute_pairwise_angles( - ... ((1.0, -1.0), (0.0, 0.0), (0.0, 0.0)), - ... )[0, 1] + >>> 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 diff --git a/test/test_dipy.py b/test/test_dipy.py index c25fedef..07b257f9 100644 --- a/test/test_dipy.py +++ b/test/test_dipy.py @@ -46,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], @@ -53,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, + ], ] ), ), @@ -125,11 +149,20 @@ ], ) def test_compute_pairwise_angles(bvecs1, bvecs2, closest_polarity, expected): - gtab1 = gradient_table([1000] * len(bvecs1), bvecs1) - gtab2 = gradient_table([1000] * len(bvecs2), bvecs2) + # DIPY requires the vectors to be normalized + _bvecs1 = bvecs1 / np.linalg.norm(bvecs1, axis=0) + gtab1 = gradient_table([1000] * _bvecs1.shape[-1], _bvecs1) + + _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) - assert (bvecs1.shape[-1], bvecs2.shape[-1]) == obtained.shape + 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) diff --git a/test/test_model.py b/test/test_model.py index 9f1fdc6e..274fabf8 100644 --- a/test/test_model.py +++ b/test/test_model.py @@ -110,7 +110,7 @@ def test_average_model(): def test_gp_model(): gp = GaussianProcessModel(kernel="default") - assert isinstance(gp, model.GaussianProcessModel) + assert isinstance(gp, model.dipy.GaussianProcessModel) X, y = make_friedman2(n_samples=500, noise=0, random_state=0) gp.fit(X, y) From 5bdff10782e123a34927b41329d3a23c2d4691ab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jon=20Haitz=20Legarreta=20Gorro=C3=B1o?= Date: Thu, 4 Jul 2024 14:06:45 -0400 Subject: [PATCH 8/9] BUG: Fix Gaussian process Fix Gaussian process: - Use a default (`test`) kernel to instantiate the Gaussian process in the test, as the spherical kernel is not yet added. - Use a regression problem to obtain 3 features (instead of 4 from the Friedman2 problem) in the test. - Normalize the gradient vectors in thest test. - Swap the signal and the gradient table when fitting the GP in the test to conform to the expected parameters. - Use the appropriate indexing and shape for the query vectors in the test. - A `dipy.core.gradients.GradientTable` does not have a `len` property so use the shape of its `gradients` attribute when checking tha the dimensionality of the data and the features match when fitting the Gaussian process. - Use named parameters when instantiating `GPFit` to avoid providing the gradient table twice. - Get the first element in the `data` vector (whose dimensionality is changed by the `GaussianProcessModel` fitting method when checking for a mask) to enable fitting. - The `model` variable in `gp_prediction` is a scikit-learn `GaussianProcessRegressor` instance do it does not have an `_gpr` attribute. --- src/eddymotion/model/dipy.py | 12 +++++++----- test/test_model.py | 16 ++++++++++------ 2 files changed, 17 insertions(+), 11 deletions(-) diff --git a/src/eddymotion/model/dipy.py b/src/eddymotion/model/dipy.py index 63e2150f..1907363c 100644 --- a/src/eddymotion/model/dipy.py +++ b/src/eddymotion/model/dipy.py @@ -64,11 +64,11 @@ def gp_prediction( # 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._gpr, "X_train_"): + 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._gpr.predict(gtab, return_std=False) + return model.predict(gtab, return_std=False) def get_kernel( @@ -175,9 +175,11 @@ def fit( data[mask[..., None]] if mask is not None else np.reshape(data, (-1, data.shape[-1])) ) - if data.shape[-1] != len(gtab): + signal_dirs = data.shape[-1] + grad_dirs = gtab.gradients.shape[0] + if signal_dirs != grad_dirs: raise ValueError( - f"Mismatched data {data.shape[-1]} and gradient table {len(gtab)} sizes." + f"Mismatched data {signal_dirs} and gradient table {grad_dirs} sizes." ) gpr = GaussianProcessRegressor( @@ -185,8 +187,8 @@ def fit( random_state=random_state, ) self._modelfit = GPFit( - gpr.fit(gtab.gradients, data), gtab=gtab, + model=gpr.fit(gtab.gradients, data[0]), mask=mask, ) return self._modelfit diff --git a/test/test_model.py b/test/test_model.py index 274fabf8..89f5235e 100644 --- a/test/test_model.py +++ b/test/test_model.py @@ -24,7 +24,8 @@ import numpy as np import pytest -from sklearn.datasets import make_friedman2 +from dipy.core.gradients import gradient_table +from sklearn.datasets import make_regression from eddymotion import model from eddymotion.data.dmri import DWI @@ -108,14 +109,17 @@ def test_average_model(): def test_gp_model(): - gp = GaussianProcessModel(kernel="default") + gp = GaussianProcessModel("test", kernel="default") assert isinstance(gp, model.dipy.GaussianProcessModel) - X, y = make_friedman2(n_samples=500, noise=0, random_state=0) - gp.fit(X, y) - X_qry = X[:2, :] - prediction, _ = gp.predict(X_qry, return_std=True) + 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],) From 5dc6c81e24f6057c18ea720849e2dd3fcc381b51 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 4 Jul 2024 22:01:41 +0200 Subject: [PATCH 9/9] Update test/test_model.py --- test/test_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_model.py b/test/test_model.py index 89f5235e..dba71522 100644 --- a/test/test_model.py +++ b/test/test_model.py @@ -109,7 +109,7 @@ def test_average_model(): def test_gp_model(): - gp = GaussianProcessModel("test", kernel="default") + gp = GaussianProcessModel("test") assert isinstance(gp, model.dipy.GaussianProcessModel)