ENH: Adds Sparse Fascicle and Gaussian Process models

[josephmje]

Addresses nipreps/nifreeze#4 and #53

[pep8speaks]

Hello @josephmje! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻

Comment last updated at 2021-12-10 20:07:16 UTC

[dPys]

In the long-term, I'd bet we get more miles out of using https://docs.pymc.io/en/v3/Gaussian_Processes.html in place of sklearn's GaussianProcessRegressor but for now this implementation should be fine!?

RE: tuning this, we could start by trying to make alpha very small like what we did with the elastic-net case, but I don't know about the kernel hyperparameters. @arokem -- would the idea here be that by wrapping this in SFM, we'd avoid the need for learning the kernel's shape? The default used for GP in sklearn is ConstantKernel(1.0, constant_value_bounds="fixed" * RBF(1.0, length_scale_bounds="fixed").

If we're sticking with this default, then I think that either way, we should remove the fixed option since "kernel hyperparameters are optimized during fitting unless the bounds are marked as fixed" and from Andersson & Sotiropoulos 2015: "The smoothness is determined from hyperparameters whose values are determined directly from the data" and "the hyperparameters will be such that the predictions are a smooth function of gradient direction and, in the case of multi-shell data, b-value." Also noteworthy was this point: "the data given by original hyperparameters is very sharp and even for a dataset with 300 points will give almost half the weight to the center-point (i.e., when predicting the signal for a diffusion direction half the information comes from the point itself)."

[dPys]

And for selecting the kernel, https://www.sciencedirect.com/science/article/pii/S1053811915006874?via%3Dihub#:~:text=Finding%20k(x,scaling%20or%20variance. they suggest using either exponential or spherical and then optimizing with maximum likelihood

[arokem]

For parallelization, we can divide the data up into chunks before calling this, and then call this separately on each chunk (on a separate thread/process)

[josephmje]

Thanks @arokem, I'm going to hand it over to @dPys to work on the chunking in a separate PR

[dPys]

@arokem and @oesteban -- Did we ever decide on whether the best approach for this would be to migrate @sebastientourbier 's asyncio parallelization routines for the DTI model to a base class for parallel fit/predict that could be imported to each of the SFM and DKI API's as well?

On a separate note-- would it be better to stick with asyncio, switch to joblib, or expose both approaches as optional?

[arokem]

I would stick with a copy of Sebastien's solution here. Let's worry about a more elegant OO design after we see that this works well here.

[oesteban]

Does this file contain b=0 volumes? I believe it should not (i.e., our dwi object doesn't filter them out)

[josephmje]

Yes it does. I think the load function filters out b values less than 50.

https://github.com/nipreps/eddymotion/blob/main/eddymotion/dmri.py#L224-L226

[josephmje]

We just realized that the load function as it currently stands affects the AverageB0 model since no b0s will be found. Will address in a separate PR.

[oesteban]

Is there an AverageB0 model? I think it's AverageDWI

[oesteban]

Meaning, it actually doesn't expect bzeros

[dPys]

LGTM!

[oesteban]

I think we should reuse Seb's solution

…

On Sat, Dec 11, 2021, 00:57 Derek Pisner ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In eddymotion/model.py <#60 (comment)>: > + 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, ...]) @arokem <https://github.com/arokem> and @oesteban <https://github.com/oesteban> -- Did we ever decide on whether the best approach for this would be to migrate @sebastientourbier <https://github.com/sebastientourbier> 's asyncio parallelization routines for the DTI model to a base class for parallel fit/predict that could be imported to each of the SFM and DKI API's as well? On a separate note-- would it be better to stick with asyncio, switch to joblib, or expose both approaches as optional? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#60 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAESDRUYYCG3SJEUTXRYDBDUQKHVDANCNFSM5JHJ4PBA> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.