From 09e5af916043a68115d69b326575779cb6ec562e Mon Sep 17 00:00:00 2001 From: Michael Joseph Date: Thu, 2 Dec 2021 09:47:31 -0500 Subject: [PATCH 1/6] wip adding gp model --- eddymotion/model.py | 74 +++++++++++++++++++++++++++++++++++++++------ 1 file changed, 65 insertions(+), 9 deletions(-) diff --git a/eddymotion/model.py b/eddymotion/model.py index b3ebe36a..6010c94b 100644 --- a/eddymotion/model.py +++ b/eddymotion/model.py @@ -25,7 +25,7 @@ def init(gtab, model="DTI", **kwargs): An array representing the gradient table in RAS+B format. model : :obj:`str` Diffusion model. - Options: ``"3DShore"``, ``"SFM"``, ``"DTI"``, ``"DKI"``, ``"S0"`` + Options: ``"3DShore"``, ``"SFM"``, ``"GP"``, ``"DTI"``, ``"DKI"``, ``"S0"`` Return ------ @@ -53,15 +53,17 @@ def init(gtab, model="DTI", **kwargs): "lambdaL": 1e-8, } - elif model.lower().startswith("sfm"): - from eddymotion.utils.model import ( - SFM4HMC as Model, - ExponentialIsotropicModel, - ) + elif model.lower() in ("sfm", "gp"): + Model = SparseFascicleModel + if model.lower() == "sfm": + param = { + "solver": "ElasticNet", + "isotropic": ExponentialIsotropicModel, + } + else: + from sklearn.gaussian_process import GaussianProcessRegressor - param = { - "isotropic": ExponentialIsotropicModel, - } + param = {"solver": GaussianProcessRegressor} elif model.lower() in ("dti", "dki"): Model = DTIModel if model.lower() == "dti" else DKIModel @@ -352,3 +354,57 @@ def _rasb2dipy(gradient): def _model_fit(model, data): return model.fit(data) + + +class SparseFascicleModel: + """ + A wrapper of :obj:`dipy.reconst.sfm.SparseFascicleModel. + """ + + __slots__ = ("_model", "_S0", "_mask", "_solver") + + def __init__(self, gtab, S0=None, mask=None, solver=None, **kwargs): + """Instantiate the wrapped model.""" + from dipy.reconst.sfm import SparseFascicleModel + from sklearn.gaussian_process import GaussianProcessRegressor + + self._S0 = None + if S0 is not None: + self._S0 = np.clip( + S0.astype("float32") / S0.max(), + a_min=1e-5, + a_max=1.0, + ) + + self._mask = mask + if mask is None and S0 is not None: + self._mask = self._S0 > np.percentile(self._S0, 35) + + if self._mask is not None: + self._S0 = self._S0[self._mask.astype(bool)] + + self._solver = solver + if solver is None: + self._solver = "ElasticNet" + + kwargs = {k: v for k, v in kwargs.items() if k in ("solver",)} + self._model = SparseFascicleModel(gtab, **kwargs) + + def fit(self, data, **kwargs): + """Clean-up permitted args and kwargs, and call model's fit.""" + self._model = self._model.fit(data[self._mask, ...]) + + def predict(self, gradient, **kwargs): + """Propagate model parameters and call predict.""" + predicted = np.squeeze( + self._model.predict( + _rasb2dipy(gradient), + S0=self._S0, + ) + ) + if predicted.ndim == 3: + return predicted + + retval = np.zeros_like(self._mask, dtype="float32") + retval[self._mask, ...] = predicted + return retval From 57db389f7f6d0da571e8e48b7c5602d77275aba4 Mon Sep 17 00:00:00 2001 From: Michael Joseph Date: Thu, 2 Dec 2021 16:03:14 -0500 Subject: [PATCH 2/6] update model --- eddymotion/model.py | 45 ++++++++++++++++++++++----------------------- 1 file changed, 22 insertions(+), 23 deletions(-) diff --git a/eddymotion/model.py b/eddymotion/model.py index 6010c94b..fb3a3f34 100644 --- a/eddymotion/model.py +++ b/eddymotion/model.py @@ -334,28 +334,6 @@ def predict(self, gradient, **kwargs): return retval -def _rasb2dipy(gradient): - gradient = np.asanyarray(gradient) - if gradient.ndim == 1: - if gradient.size != 4: - raise ValueError("Missing gradient information.") - gradient = gradient[..., np.newaxis] - - if gradient.shape[0] != 4: - gradient = gradient.T - elif gradient.shape == (4, 4): - print("Warning: make sure gradient information is not transposed!") - - with warnings.catch_warnings(): - warnings.filterwarnings("ignore", category=UserWarning) - retval = gradient_table(gradient[3, :], gradient[:3, :].T) - return retval - - -def _model_fit(model, data): - return model.fit(data) - - class SparseFascicleModel: """ A wrapper of :obj:`dipy.reconst.sfm.SparseFascicleModel. @@ -366,7 +344,6 @@ class SparseFascicleModel: def __init__(self, gtab, S0=None, mask=None, solver=None, **kwargs): """Instantiate the wrapped model.""" from dipy.reconst.sfm import SparseFascicleModel - from sklearn.gaussian_process import GaussianProcessRegressor self._S0 = None if S0 is not None: @@ -408,3 +385,25 @@ def predict(self, gradient, **kwargs): retval = np.zeros_like(self._mask, dtype="float32") retval[self._mask, ...] = predicted return retval + + +def _rasb2dipy(gradient): + gradient = np.asanyarray(gradient) + if gradient.ndim == 1: + if gradient.size != 4: + raise ValueError("Missing gradient information.") + gradient = gradient[..., np.newaxis] + + if gradient.shape[0] != 4: + gradient = gradient.T + elif gradient.shape == (4, 4): + print("Warning: make sure gradient information is not transposed!") + + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", category=UserWarning) + retval = gradient_table(gradient[3, :], gradient[:3, :].T) + return retval + + +def _model_fit(model, data): + return model.fit(data) From fc771f58da00c28f13ad8ff6382048378fe8bb68 Mon Sep 17 00:00:00 2001 From: Michael Joseph Date: Fri, 10 Dec 2021 09:03:40 -0500 Subject: [PATCH 3/6] Apply suggestions from code review Co-authored-by: Oscar Esteban --- eddymotion/model.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/eddymotion/model.py b/eddymotion/model.py index fb3a3f34..27a9e310 100644 --- a/eddymotion/model.py +++ b/eddymotion/model.py @@ -55,12 +55,11 @@ def init(gtab, model="DTI", **kwargs): elif model.lower() in ("sfm", "gp"): Model = SparseFascicleModel - if model.lower() == "sfm": - param = { - "solver": "ElasticNet", - "isotropic": ExponentialIsotropicModel, - } - else: + param = { + "solver": "ElasticNet", + "isotropic": ExponentialIsotropicModel, + } + if model.lower() == "gp": from sklearn.gaussian_process import GaussianProcessRegressor param = {"solver": GaussianProcessRegressor} From b4624da3e52174fe876cdce5e3eaf1f055c170e9 Mon Sep 17 00:00:00 2001 From: Michael Joseph Date: Fri, 10 Dec 2021 10:07:27 -0500 Subject: [PATCH 4/6] update requirements with missing packages --- setup.cfg | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index 3727bba2..96e9413c 100644 --- a/setup.cfg +++ b/setup.cfg @@ -22,8 +22,11 @@ url = https://github.com/nipreps/EddyMotionCorrection python_requires = >=3.7 install_requires = dipy>=1.3.0 - scikit-image>=0.14.2 + nipype>= 1.5.1, < 2.0 + nitransforms>=21.0.0 nest-asyncio>=1.5.1 + scikit-image>=0.14.2 + scikit-learn>=1.0.1 test_requires = codecov coverage From 5f6dd9ddad8ac17d0f06519fbe508a332b7498fe Mon Sep 17 00:00:00 2001 From: Michael Joseph Date: Fri, 10 Dec 2021 11:05:51 -0500 Subject: [PATCH 5/6] add ExponentialIsotropicModel if multi-shell --- eddymotion/model.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/eddymotion/model.py b/eddymotion/model.py index 27a9e310..f990ef61 100644 --- a/eddymotion/model.py +++ b/eddymotion/model.py @@ -55,10 +55,14 @@ def init(gtab, model="DTI", **kwargs): elif model.lower() in ("sfm", "gp"): Model = SparseFascicleModel - param = { - "solver": "ElasticNet", - "isotropic": ExponentialIsotropicModel, - } + param = {"solver": "ElasticNet"} + + from dipy.core.gradients import check_multi_b + multi_b = check_multi_b(gtab, 2, non_zero=False) + if multi_b: + from dipy.reconst.sfm import ExponentialIsotropicModel + param.update({"isotropic": ExponentialIsotropicModel}) + if model.lower() == "gp": from sklearn.gaussian_process import GaussianProcessRegressor From b232df697520ff71e2c526c2889650b20b446da4 Mon Sep 17 00:00:00 2001 From: Michael Joseph Date: Fri, 10 Dec 2021 15:07:09 -0500 Subject: [PATCH 6/6] swap order of assigning isotropic key --- eddymotion/model.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/eddymotion/model.py b/eddymotion/model.py index f990ef61..ed0725b2 100644 --- a/eddymotion/model.py +++ b/eddymotion/model.py @@ -6,7 +6,7 @@ import nest_asyncio import numpy as np -from dipy.core.gradients import gradient_table +from dipy.core.gradients import check_multi_b, gradient_table nest_asyncio.apply() @@ -57,17 +57,15 @@ def init(gtab, model="DTI", **kwargs): Model = SparseFascicleModel param = {"solver": "ElasticNet"} - from dipy.core.gradients import check_multi_b + if model.lower() == "gp": + from sklearn.gaussian_process import GaussianProcessRegressor + param = {"solver": GaussianProcessRegressor} + multi_b = check_multi_b(gtab, 2, non_zero=False) if multi_b: from dipy.reconst.sfm import ExponentialIsotropicModel param.update({"isotropic": ExponentialIsotropicModel}) - if model.lower() == "gp": - from sklearn.gaussian_process import GaussianProcessRegressor - - param = {"solver": GaussianProcessRegressor} - elif model.lower() in ("dti", "dki"): Model = DTIModel if model.lower() == "dti" else DKIModel