Return S0 value for DTI fits #1162
Conversation
Added tests for fast linear fitting in test_ivim.py Added a fitting comparison function that uses AIC
Just getting it up to date
Also added bounds checking and clamping options. Also made it so the relative likelihood can be positive or negative.
think it's working correctly. Fixed errors in the bound checking and clipping.
Changed test accuracy parameter.
Added an AIC weight calculation function.
Added tests for the aic functions.
Changed the tests to reflect this.
Added the requisite tests.
Fixed a bug that caused an error when asking for the S0 values when they weren't calculated.
Removed the error for the above methods preventing S0 fitting in those cases.
…ted. Fixed the S0 prediction test.
| S0params[i:i + step] = fit_output[1] | ||
| else: | ||
| dtiparams[i:i + step] = fit_tensor(design_matrix, | ||
| data[i:i + step], *args, **kwargs) |
There was a problem hiding this comment.
PEP8: align data[i:i + step]
There was a problem hiding this comment.
Also, ensure that you don't have a line with more than 79 characters
|
That sounds reasonable to me. I made the changes so S0_hat is the model fit. I think it's working and doing everything it should. I suspect I'm missing some test cases and documentation so I'll try to check that when I can but comments are welcome before then. |
Changed the internal S0_hat variable to model_S0 for consistency. Added S0_hat in TensorFit __getitem__
|
I updated the tests and documentation. Let me know what you think. I'm also interested if you have thoughts about speed and memory testing beyond the basic script I'm using locally. |
There was a problem hiding this comment.
Nice work -- thanks for implementing this! I wonder whether some of the logic could be simplified a bit (see comments inline). Regarding profiling: you can use this script: https://github.com/nipy/dipy/blob/master/scratch/profile_dti.py, together with the line-profiler: https://github.com/rkern/line_profiler to measure timing, and profile this line-by-line. At the very least using %timeit to compare performance with and without this would be good, just to see that we are not slowing this down for people who don't care about estimating S0.
One more request I have is to test this with different configurations of b=0 measurements: with one measurement at the beginning, and at the end of the acquisition, with several b=0 measurements, either bunched together, or interspersed throughout. I want to be sure that this is not somehow making assumptions about the b=0 measurements that just happen to be true for the test conditions.
| params_in_mask = self.fit_method(self.design_matrix, data_in_mask, | ||
| *self.args, **self.kwargs) | ||
|
|
||
| if self.return_S0_hat: |
There was a problem hiding this comment.
I don't think you need this if... else construct. If it's false, it will just do the right thing (assuming it's implemented correctly in the fit_method).
There was a problem hiding this comment.
I agree on the input side that's fine but the output of fit_method is a little more tricky because they either return 1 or 2 arguments. So I think the alternative below still requires an if statement but is perhaps more clean?
params_in_mask = self.fit_method(
self.design_matrix,
data_in_mask,
return_S0_hat=self.return_S0_hat,
*self.args,
**self.kwargs)
if self.return_S0_hat:
params_in_mask, model_S0 = params_in_maskThere was a problem hiding this comment.
Yes - that's cleaner. In general, I prefer it when there is only one call to the function. That way, if the API changes for reasons unrelated to S0_hat, we only have to change the call in one place, instead of two places.
|
|
||
| if mask is None: | ||
| out_shape = data.shape[:-1] + (-1, ) | ||
| dti_params = params_in_mask.reshape(out_shape) | ||
| if self.return_S0_hat: |
There was a problem hiding this comment.
I think the logic here could be simplified a bit by setting model_S0 to None, changing it if required, and otherwise passing it on as None to TensorFit below.
There was a problem hiding this comment.
Could you explain how to simplify the if statements in the mask logic? I don't see how that would work. If model_S0 is None then it would fail trying to reshape or it would return a large None array.
I do agree that the ifs with TensorFit can be removed though to leave return TensorFit(self, dti_params, model_S0=S0_params) as long as S0_params is None.
There was a problem hiding this comment.
Yes - I would make S0_params an optional input to fit with a default of None. Does that make sense? But you're right that I don't see a way to simplify more.
| N = model_params.ndim | ||
| if type(index) is not tuple: | ||
| index = (index,) | ||
| elif len(index) >= model_params.ndim: | ||
| raise IndexError("IndexError: invalid index") | ||
| index = index + (slice(None),) * (N - len(index)) | ||
| return type(self)(self.model, model_params[index]) | ||
| if model_S0 is None: |
There was a problem hiding this comment.
And same here. If None was passed through, you can pass probably pass it on to the call to type(self)(self.model, ...
There was a problem hiding this comment.
But model_S0=model_S0[index_S0] would fail if model_S0 is None. Or are you suggesting just pass it along if it's None anyway on line 848?
There was a problem hiding this comment.
Something like:
if model_S0 is not None:
model_S0 = model_S0[index[:-1]]
return type(self)(self.model, model_params[index], model_S0=model_S0)
| @@ -1162,6 +1196,10 @@ def predict(self, gtab, S0=1., step=None): | |||
| which a signal is to be predicted and $b$ is the b value provided in | |||
| the GradientTable input for that direction | |||
| """ | |||
| if S0 is None: | |||
There was a problem hiding this comment.
This should be documented in the docstring. In particular, now that the "1" is now no longer part of the function signature, it is worth mentioning that value specifically in the docstring for this function.
| @@ -1242,22 +1283,41 @@ def wrapped_fit_tensor(design_matrix, data, step=step, | |||
| size = np.prod(shape) | |||
| step = int(step) or size | |||
| if step >= size: | |||
| return fit_tensor(design_matrix, data, *args, **kwargs) | |||
| if return_S0_hat: | |||
There was a problem hiding this comment.
Here as well, consider whether the logic could be simplified by passing along the default value in the case where return_S0_hat is set to False.
There was a problem hiding this comment.
You're right, this can be simplified. Though not anywhere else in this function as far as I can see.
| @@ -631,6 +657,11 @@ def test_restore(): | |||
| tensor_model = dti.TensorModel(gtab, fit_method='restore', sigma=0.0001) | |||
| tensor_model.fit(Y.copy()) | |||
|
|
|||
| # Test return_S0_hat | |||
| tensor_model = dti.TensorModel(gtab, fit_method='restore', sigma=0.0001, return_S0_hat=True) | |||
|
I'll give profiling a shot. I'll try to do this PR with and without returning S0 and master. I'll also see if I can track memory in a meaningful way. Regarding testing different structured inputs, is there a dataset that's already in dipy you recommend? Do you think it would be good to test on data without any b=0 at all? |
|
Here are some times, timed with the version that's up right now. It's not taking much more time at all, but the memory size does increase because we need to hold on to the S0 data rather than throw it away. Master for comparison because it has no S0 logic at all: Current version, S0 OFF: Current version, S0 ON: |
|
@arokem your comments sound good to me so I made and pushed the changes. Also, do you want more details about timing and profiling or is what I have above using timeit satisfactory? |
|
Code looks good to me. Regarding profiling: what data did you use to run the profiling? Assuming this is the difference on reasonably large whole-brain data, it seems to be that there wouldn't be harm in adding this code -- the difference seems minor. Any thoughts from anyone else? |
|
Good to hear. I used some 3T human brain data I have here that has multiple b-values so the S0 fitting is meaningful. It's 256x256x64 and 3 directions with 7 b-values. I did mask it but it's still fitting 905,768 voxels. |
|
A somewhat related question: is Dipy participating in Google Summer of Code this year? I'm thinking about this PR and integrating it into IVIM and other IVIM improvements that could work as a project for GSoC. I think the first deadline is Feb 9th so it's already looming on the horizon. |
|
OK - I am happy with this. If no one else has any comments here, I will merge this in a couple of days. Re: GSoC -- I will not be able to be involved this year. I am going to be on family leave during the summer 🍼 😄 ! |
|
Congrats on the 👶 ! That's a legit reason to skip it. We're expecting our first in February so I was figuring I'd at least be partly back in action by April in time for GSoC. But maybe less is more this year 😄 |
|
I've just did a last look to the code - for me this is ready to go! Thanks @etpeterson for the really nice work! |
Return S0 value for DTI fits
The goal of this PR is to optionally return the S0 (b=0) value for DTI fits without negatively impacting performance or memory unless the S0 is requested. It captures the fitted S0 values generated in the DTI fitting routines and allows the user to access them.
It's similar in intent to PR 1036 but rather than calculating the S0 values from the fits, it captures the values directly in each fit method.
My local tests look like the speed is only slightly slower but keeping the S0 values does increase memory usage because it's another image being saved. I'm timing it using time for the whole fit and finding the size by pickling the fit object so if there's a better way to do that let me know.
I think it's missing a few tests and probably some documentation but I think the major pieces are there.