ENH: Nipype-less implementation (refactor) #16
Conversation
There was a problem hiding this comment.
Looks good! One general point is that I think that this probably ends up looking more like sklearn unsupervised estimators. For example, note that the predict method takes no input. This means that it is probably more akin to the sklearn transform methods in unsupervised estimators. And we can add a fit_transform that does both. Other small comments spread throughout.
emc/dmri.py
Outdated
| class DWI: | ||
| """Data representation structure for dMRI data.""" | ||
|
|
||
| dataobj = attr.ib(default=None) |
There was a problem hiding this comment.
Does this name follow a specific convention? Maybe better to call it dwdata to be more explicit that these are the diffusion weighted images?
There was a problem hiding this comment.
nibabel gives access to the numpy array through dataobj, but it is not required, honestly.
emc/estimator.py
Outdated
| gradients = attr.ib(default=None) | ||
| """A 2D numpy array of the gradient table in RAS+B format.""" | ||
| bzero = attr.ib(default=None) | ||
| """A *b=0* reference map, preferably obtained by some smart averaging.""" |
There was a problem hiding this comment.
Is this a path or an array? I would advocate for either.
|
I hadn't assigned review because it felt premature - but this feedback is super-welcome. Thanks for your time, I will move to the better API tomorrow. |
|
@arokem I'm getting this: >>> EddyMotionEstimator.fit(data)
Optimizing level 0 [max iter: 100]
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-35-e34da2fc4915> in <module>
----> 1 EddyMotionEstimator.fit(data)
~/workspace/EddyMotionCorrection/emc/estimator.py in fit(dwdata, n_iter, align_kwargs, model, seed, **kwargs)
99
100 init_affine = dwdata.em_affines[i] if dwdata.em_affines else np.eye(4)
--> 101 xform = registration.optimize(
102 predicted, # fixed
103 data_test[0], # moving
~/.anaconda/lib/python3.8/site-packages/dipy/align/imaffine.py in optimize(***failed resolving arguments***)
1099 self.options['maxiter'] = max_iter
1100
-> 1101 opt = Optimizer(self.metric.distance_and_gradient,
1102 self.params0,
1103 method=self.method, jac=True,
~/.anaconda/lib/python3.8/site-packages/dipy/core/optimize.py in __init__(self, fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options, evolution)
127 else:
128
--> 129 res = minimize(fun, x0, args, method, jac, hess, hessp, bounds,
130 constraints, tol, callback, options)
131
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
615 **options)
616 elif meth == 'l-bfgs-b':
--> 617 return _minimize_lbfgsb(fun, x0, args, jac, bounds,
618 callback=callback, **options)
619 elif meth == 'tnc':
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/lbfgsb.py in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, finite_diff_rel_step, **unknown_options)
304 iprint = disp
305
--> 306 sf = _prepare_scalar_function(fun, x0, jac=jac, args=args, epsilon=eps,
307 bounds=new_bounds,
308 finite_diff_rel_step=finite_diff_rel_step)
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/optimize.py in _prepare_scalar_function(fun, x0, jac, args, bounds, epsilon, finite_diff_rel_step, hess)
259 # ScalarFunction caches. Reuse of fun(x) during grad
260 # calculation reduces overall function evaluations.
--> 261 sf = ScalarFunction(fun, x0, args, grad, hess,
262 finite_diff_rel_step, bounds, epsilon=epsilon)
263
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py in __init__(self, fun, x0, args, grad, hess, finite_diff_rel_step, finite_diff_bounds, epsilon)
74
75 self._update_fun_impl = update_fun
---> 76 self._update_fun()
77
78 # Gradient evaluation
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py in _update_fun(self)
164 def _update_fun(self):
165 if not self.f_updated:
--> 166 self._update_fun_impl()
167 self.f_updated = True
168
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py in update_fun()
71
72 def update_fun():
---> 73 self.f = fun_wrapped(self.x)
74
75 self._update_fun_impl = update_fun
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py in fun_wrapped(x)
68 def fun_wrapped(x):
69 self.nfev += 1
---> 70 return fun(x, *args)
71
72 def update_fun():
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/optimize.py in __call__(self, x, *args)
72 def __call__(self, x, *args):
73 """ returns the the function value """
---> 74 self._compute_if_needed(x, *args)
75 return self._value
76
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/optimize.py in _compute_if_needed(self, x, *args)
66 if not np.all(x == self.x) or self._value is None or self.jac is None:
67 self.x = np.asarray(x).copy()
---> 68 fg = self.fun(x, *args)
69 self.jac = fg[1]
70 self._value = fg[0]
~/.anaconda/lib/python3.8/site-packages/dipy/align/imaffine.py in distance_and_gradient(self, params)
783 """
784 try:
--> 785 self._update_mutual_information(params, True)
786 except (AffineInversionError, AffineInvalidValuesError):
787 return np.inf, 0 * self.metric_grad
~/.anaconda/lib/python3.8/site-packages/dipy/align/imaffine.py in _update_mutual_information(self, params, update_gradient)
670
671 # Update the histogram with the current joint intensities
--> 672 static_values, moving_values = self._update_histogram()
673
674 H = self.histogram # Shortcut to `self.histogram`
~/.anaconda/lib/python3.8/site-packages/dipy/align/imaffine.py in _update_histogram(self)
631 pts = sp_to_moving.dot(self.samples.T).T # Points on moving grid
632 pts = pts[..., :self.dim]
--> 633 self.moving_vals, inside = self.interp_method(self.moving, pts)
634 self.moving_vals = np.array(self.moving_vals)
635 static_values = self.static_vals
dipy/core/interpolation.pyx in dipy.core.interpolation.__pyx_fused_cpdef()
TypeError: No matching signature found |
|
Okay, nevermind for now @arokem. I have moved on to
|
|
Okay, I think we can get this merged (when reviews show a thumbs up) and then see where we go from here (basically, use a non-dummy prediction model) |
|
This does significantly change the objectives of the software, which was a pure python implementation of these ideas. Regarding API/usage: I think we can sort that out. Making it easier for newbies should be a matter of writing clearer APIs. That should also help you as a developer, I believe. Regarding finesse: if these features are important for performance, we can implement this in Python. Regarding division of responsibility: I am not sure I follow. What's the issue with relying on DIPY? |
|
I guess one way we could go about this is to merge this now, and then it would be up to me to refactor this out later on. I'd be OK with that. |
|
One issue for the tutorial (and in general) : how do you expect users to install the software? A |
|
Yes, I agree this raises a whole lot of questions and I also want to see this being python-only. Conda could be a good trade-off for installation until that happens. For the tutorial, a docker image should do. |
|
Re: division of responsibilities If we feel there are some pieces of code that would make sense to upstream into DIPY, it will be easier to do if we focus all the attention (about DIPY) on the models. When that bit is rock-solid, we can then proceed with the registration part. Also, getting good results with ants can be the best way to set a reference quality for DIPY's implementation. |
|
I would like to merge this asap - may I ask for a final round of reviews? @arokem @dPys @sebastientourbier |
There was a problem hiding this comment.
Such incredible progress on this! I think users are really going to love the sklearn-style fit, transform, predict methods. I left a few comments and a suggested enhancement, but otherwise just very excited to see this tool maturing, especially with the "Oscar-touch" :-)
|
|
||
| with h5py.File(filename, "w") as out_file: | ||
| out_file.attrs["Format"] = "EMC/DWI" | ||
| out_file.attrs["Version"] = np.uint16(1) |
There was a problem hiding this comment.
Can we add more metadata to this? Now that we're moving towards HDF5 and away from MapNodes, It would be insanely helpful to tag each image file that gets written out with its provenance (is this all volumes?, DWI with or without its B0 volume(s)?, masked or unmasked? raw, transformed or predicted? etc. If it's transformed, is it an eddy-corrected transformation or an eddy->predicted transformation. If it was predicted, what base model was used to predict it? Coupled with the rasb gradient information as well, it might be very helpful to have this class automatically track the indices of individual volumes of the images themselves. WDYT?
There was a problem hiding this comment.
Actually, this file contains the original dataset. I'm still unsure of the I/O because the current implementation allows you to do some funky stuff (e.g., if you change the _filepath attribute, the behavior is easily problematic because you could potentially store a realigned dataset and that's not what the estimator is expecting).
But yes, happy to grow from here.
|
|
||
| param = { | ||
| "isotropic": ExponentialIsotropicModel, | ||
| } |
There was a problem hiding this comment.
Can we go ahead and add alpha + lambda to this list?
There was a problem hiding this comment.
Let's get this merged and sure, on top of this PR we can definitely do that.
| index_order = np.arange(len(dwdata)) | ||
| np.random.shuffle(index_order) | ||
| with tqdm(total=len(index_order), unit="dwi") as pbar: | ||
| for i in index_order: |
There was a problem hiding this comment.
What about starting out with just joblib?
There was a problem hiding this comment.
That would parallelize whole orientations. @sebastientourbier is preparing the parallelization of blocks within any given orientation. The idea here is to be able to benefit from previous alignments in the latter alignments.
| # read original DWI data & b-vector | ||
| with h5py.File(self._filepath, "r") as in_file: | ||
| root = in_file["/0"] | ||
| dwframe = np.asanyarray(root["dataobj"][..., index]) |
There was a problem hiding this comment.
curious -- why dwframe and bframe vs. just test_data and test_gradients?
There was a problem hiding this comment.
test_data and test_gradients imply there is some cross-validation going on. Because this is just a data representation object, it seemed more adequate not to refer to the framework you are working on.
|
|
||
| train_data = self.dataobj[..., ~mask] | ||
| train_gradients = self.gradients[..., ~mask] | ||
|
|
There was a problem hiding this comment.
So I think it would still be good to have some mechanism for further filtering which gradient samples are included for model training and testing upon each iteration. On the one hand, there are those samples that are too close in proximity along the sphere to those which we are trying to predict and therefore risk unfavorably biasing the predictions. But even if that's still a matter of debate, there may simply be outlier volumes or gradients that we wish to manually exclude. For that reason alone, I think it would still be beneficial to assume some additional level of control over which samples get used, no?
On that note, I do think that a threshold for filtering out training samples based on q-space proximity should itself be adaptive since I'd assume that the extent of this bias will vary as a function of model choices and sampling scheme (e.g. we should also think about what might happen here in the half-sphere case and how best to accommodate...).
But for now, adding some functionality that does something like the following might be a good starting point:
def train_gradient_filter(
test_gradients, train_gradients, cutoff=2.0
):
train_bvec = np.asanyarray(train_gradients[:, 0:3])
train_bval = np.asanyarray(train_gradients[:, 3])
test_bvec = np.asanyarray(test_gradients[:, 0:3])
test_bval = np.asanyarray(test_gradients[:, 3])
min_bval = min(min(train_bval), test_bval)
all_qvals = np.sqrt(train_bval - min_bval)
prediction_qval = np.sqrt(test_bval - min_bval)
# Convert q values to percent of maximum qval
max_qval = max(max(all_qvals), prediction_qval)
all_qvals_scaled = all_qvals / max_qval * 100
scaled_qvecs = train_bvec * all_qvals_scaled[:, np.newaxis]
scaled_prediction_qvec = test_bvec * (prediction_qval / max_qval * 100)
# Calculate the distance between the training qvecs and the testing qvec
gradient_mask = (
np.linalg.norm(scaled_qvecs - scaled_prediction_qvec, axis=1
) > cutoff
) * (np.linalg.norm(scaled_qvecs + scaled_prediction_qvec, axis=1
) > cutoff)
return gradient_mask
gradient_mask = train_gradient_filter(test_gradients, train_gradients)
train_data = self.dataobj[..., ~gradient_mask]
train_gradients = self.gradients[..., ~gradient_mask]
And then maybe also adding an exclusion_indices argument as well for manually dropping training samples if needed.
WDYT @oesteban ?
There was a problem hiding this comment.
therefore risk unfavorably biasing the predictions.
It is still unclear to me why the bias would be unfavorable in this case
there may simply be outlier volumes or gradients that we wish to manually exclude.
This is a chicken-and-egg problem. How do you decide what volumes should be excluded? Perhaps they'd be acceptable after head motion? I prefer to make the implementation as agnostic as possible of the rest of the pipeline. This is not to say that we should not look into this issue, but my opinion is that we can only investigate this outlier filtering after we have successfully tried the code on a lot of data. We need to measure how much it improves w.r.t. a baseline. This PR tries to set such a baseline.
Instead of thinking about solutions to hypothetical problems, let's hypothesize about the most relevant problems we expect first, then plan for a test that would evaluate solutions to the problem, and finally once we have both hypothesis and test, make decisions to improve performance.
Finally, and in particular about the code you propose:
- I would not implement this in the DWI object because this is not partitioning the dataset, it is partioning the partition with some additional criteria. IMHO, the responsibility of this falls more on the estimator object side.
- Let's keep this code around and test it when we have a functional prototype (i.e., a baseline as I called it before). This first implementation only looked for setting up the foundations.
There was a problem hiding this comment.
For instance, one problem I do see coming at us:
Say you have 12 b=1000, 12 b=2000 and 12 b=3000. Eddy currents for the two outer shells will be much larger than those of the b=1000. Because the dataset is already pretty biased towards the outer shells, when aligning one orientation of the b=1000 shell we might actually "overcorrect" for eddy. So in this case, you might want to completely dismiss the largest b (instead of close orientations) when "predicting" the lowest b. WDYT @arokem ?
Co-authored-by: Ariel Rokem <arokem@gmail.com>
I'm making progress on the integration with DIPY's regitration. However,
I'm hitting the following issue:
```Python
>>> EddyMotionEstimator.fit(data)
Optimizing level 0 [max iter: 100]
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-35-e34da2fc4915> in <module>
----> 1 EddyMotionEstimator.fit(data)
~/workspace/EddyMotionCorrection/emc/estimator.py in fit(dwdata, n_iter, align_kwargs, model, seed, **kwargs)
99
100 init_affine = dwdata.em_affines[i] if dwdata.em_affines else np.eye(4)
--> 101 xform = registration.optimize(
102 predicted, # fixed
103 data_test[0], # moving
~/.anaconda/lib/python3.8/site-packages/dipy/align/imaffine.py in optimize(***failed resolving arguments***)
1099 self.options['maxiter'] = max_iter
1100
-> 1101 opt = Optimizer(self.metric.distance_and_gradient,
1102 self.params0,
1103 method=self.method, jac=True,
~/.anaconda/lib/python3.8/site-packages/dipy/core/optimize.py in __init__(self, fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options, evolution)
127 else:
128
--> 129 res = minimize(fun, x0, args, method, jac, hess, hessp, bounds,
130 constraints, tol, callback, options)
131
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
615 **options)
616 elif meth == 'l-bfgs-b':
--> 617 return _minimize_lbfgsb(fun, x0, args, jac, bounds,
618 callback=callback, **options)
619 elif meth == 'tnc':
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/lbfgsb.py in _minimize_lbfgsb(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, finite_diff_rel_step, **unknown_options)
304 iprint = disp
305
--> 306 sf = _prepare_scalar_function(fun, x0, jac=jac, args=args, epsilon=eps,
307 bounds=new_bounds,
308 finite_diff_rel_step=finite_diff_rel_step)
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/optimize.py in _prepare_scalar_function(fun, x0, jac, args, bounds, epsilon, finite_diff_rel_step, hess)
259 # ScalarFunction caches. Reuse of fun(x) during grad
260 # calculation reduces overall function evaluations.
--> 261 sf = ScalarFunction(fun, x0, args, grad, hess,
262 finite_diff_rel_step, bounds, epsilon=epsilon)
263
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py in __init__(self, fun, x0, args, grad, hess, finite_diff_rel_step, finite_diff_bounds, epsilon)
74
75 self._update_fun_impl = update_fun
---> 76 self._update_fun()
77
78 # Gradient evaluation
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py in _update_fun(self)
164 def _update_fun(self):
165 if not self.f_updated:
--> 166 self._update_fun_impl()
167 self.f_updated = True
168
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py in update_fun()
71
72 def update_fun():
---> 73 self.f = fun_wrapped(self.x)
74
75 self._update_fun_impl = update_fun
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py in fun_wrapped(x)
68 def fun_wrapped(x):
69 self.nfev += 1
---> 70 return fun(x, *args)
71
72 def update_fun():
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/optimize.py in __call__(self, x, *args)
72 def __call__(self, x, *args):
73 """ returns the the function value """
---> 74 self._compute_if_needed(x, *args)
75 return self._value
76
~/.anaconda/lib/python3.8/site-packages/scipy/optimize/optimize.py in _compute_if_needed(self, x, *args)
66 if not np.all(x == self.x) or self._value is None or self.jac is None:
67 self.x = np.asarray(x).copy()
---> 68 fg = self.fun(x, *args)
69 self.jac = fg[1]
70 self._value = fg[0]
~/.anaconda/lib/python3.8/site-packages/dipy/align/imaffine.py in distance_and_gradient(self, params)
783 """
784 try:
--> 785 self._update_mutual_information(params, True)
786 except (AffineInversionError, AffineInvalidValuesError):
787 return np.inf, 0 * self.metric_grad
~/.anaconda/lib/python3.8/site-packages/dipy/align/imaffine.py in _update_mutual_information(self, params, update_gradient)
670
671 # Update the histogram with the current joint intensities
--> 672 static_values, moving_values = self._update_histogram()
673
674 H = self.histogram # Shortcut to `self.histogram`
~/.anaconda/lib/python3.8/site-packages/dipy/align/imaffine.py in _update_histogram(self)
631 pts = sp_to_moving.dot(self.samples.T).T # Points on moving grid
632 pts = pts[..., :self.dim]
--> 633 self.moving_vals, inside = self.interp_method(self.moving, pts)
634 self.moving_vals = np.array(self.moving_vals)
635 static_values = self.static_vals
dipy/core/interpolation.pyx in dipy.core.interpolation.__pyx_fused_cpdef()
TypeError: No matching signature found
```
|
Okay, we will need to push this faster. Happy to discuss all the lines above in future PRs building on top of this one. Thanks for the feedback! |
ENH: Nipype-less implementation (refactor)
ENH: Nipype-less implementation (refactor)
ENH: Nipype-less implementation (refactor) Former-commit-id: 57c5189
Setting some foundations for the code:
UPDATE:
This is how it looks right now: