Xval #335

arokem · 2014-03-11T23:37:14Z

This revives #235. Now that all the other elements are in place, I think that this is ready for review and in principle ready to be merged upstream (once it's been reviewed and fixed up, of course!). I am excited about this PR - evaluation of DWI models is something that I have been thinking a lot about in the past couple of years, since I started working on DWI.

This PR also demonstrates the power and elegance of the uniform Model/Fit API we have in place. GGT!

For now, this PR only implements cross-validation for CSD and for the tensor model, but future work could extend this to other models, by implementing their predict methods.

Also : a lot of PEP8 fixes in the CSD module and elsewhere :-)

Garyfallidis · 2014-03-17T00:32:00Z

Wow! This is cool. Have you tested what happens when a model does not provide a predict method?
And was it really necessary to use the CSDFit rather than SHFit for prediction in csdeconv? Would it be easy to have the same prediction for all SH-like models? Just saying... :)

arokem · 2014-03-17T16:21:55Z

I think that if it doesn't have a predict method, it should throw an error. I will implement that. Concerning the SHFit - I initially tried that, but couldn't make it work - but I agree that it would be better - I will give it another try.

arokem · 2014-03-17T18:57:42Z

OK - these commits should answer @Garyfallidis questions. The issue with other SHM-based models is that they don't have a response function defined, so we need to fall back on a default response function (here it's a tensor with evals [0.0015, 0.0003, 0.0003]). Does this look OK to you?

MrBago · 2014-03-18T02:29:13Z

@arokem, this stuff looks awesome. I don't think generalizing the idea of a response function to all odf models is conceptually sound, but think we should try and implement predict methods for as many models/fits as possible.

I think for all our predict methods we should include a closed loop test, IE fit(predict(P)) == P, where P is some set of parameters to the of the model. I'd like to hear if you guys think there is a good reason some models cannot pass this closed loop test.

I think, both from a design point of view and a conceptual point of view, predict should defer to the model. Somehow it makes sense that the models know the best way to get from a fit or a set of parameters back to the signal. Also this will allow us to use the same fit type for multiple models like we're currently doing with the SphHarmFit. That's not to say we can't have a fit.predict method, I really like the syntax:

model = CsaOdfMode(...)
fit = model.fit(data)
pred_signal = fit.predict()

arokem · 2014-03-18T18:06:07Z

Thanks for the feedback @MrBago!

On Mon, Mar 17, 2014 at 7:29 PM, MrBago notifications@github.com wrote:

@arokem https://github.com/arokem, this stuff looks awesome. I don't
think generalizing the idea of a response function to all odf models is
conceptually sound, but think we should try and implement predict methods
for as many models/fits as possible.

Agreed! I think that this is only the first step.

Does it make sense to have the response function in the CSA model? Right
now, the implementation of the "ODF convolved with response function" idea
is at the level of the to the base SHM model Fit object. Does that make
sense?

I think for all our predict methods we should include a closed loop test,
IE fit(predict(P)) == P, where P is some set of parameters to the of the
model. I'd like to hear if you guys think there is a good reason some
models cannot pass this closed loop test.

One reason is regularization. Some models (CSD is one example) would not do
that. As you can see in the test I implemented, CSD doesn't predict
perfectly even in a simulation when the exact response function is known
and when there is no noise. I believe this is because of the regularization
implemented in this model.

I think, both from a design point of view and a conceptual point of view,
predict should defer to the model. Somehow it makes sense that the models
know the best way to get from a fit or a set of parameters back to the
signal. Also this will allow us to use the same fit type for multiple
models like we're currently doing with the SphHarmFit. That's not to say
we can't have a fit.predict method, I really like the syntax:

model = CsaOdfMode(...)
fit = model.fit(data)
pred_signal = fit.predict()

Yeah - that's the uniform API that makes the cross-validation part of this
PR so general, that it could (in principle) be applied to any model!

But you think that this method should call out to a method of the Model
object? I am not sure that makes sense. A prediction can only really be
made from a fit, because it requires model parameters. So - a model by
itself would know something about the conceptual structure (which is why
you might need to query self.model attributes within fit.predict), but
it doesn't know how to predict, unless you provide it with parameters.
Alternatively, we could implement an alternative model.predict that takes
parameters as input and does the prediction. Is that what you were
thinking? In a way, it's similar to what we do in the sims module.

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/335#issuecomment-37893885
.

MrBago · 2014-03-20T23:26:17Z

dipy/reconst/odf.py

@@ -20,6 +22,9 @@ def fit(self, data):

 class OdfFit(ReconstFit):

+    def __init__(self, model, data):
+        ReconstFit.__init__(self, model, data)


Shouldn't OdfFit inherit init from ReconstFit automatically? Is there a reason you're doing this kind of explicit inheritance?

MrBago · 2014-03-21T01:00:24Z

Does it make sense to have the response function in the CSA model? Right
now, the implementation of the "ODF convolved with response function" idea
is at the level of the to the base SHM model Fit object. Does that make
sense?

We're abusing the ODF concept a little, in principle the csd model produces
an FOD not an ODF,
but for API consistency we're calling it an ODF. It doesn't may any
theoretical sense to convolve
and ODF with a response function, but I don't know if anyone has tried it.
Maybe it does something nice.

One reason is regularization. Some models (CSD is one example) would not
do
that. As you can see in the test I implemented, CSD doesn't predict
perfectly even in a simulation when the exact response function is known
and when there is no noise. I believe this is because of the
regularization
implemented in this model.

The CSD model should pass the closed loop test for well behaved signals,
meaning fit(predict(P)) == P if all(fod(P) > 0) (does that make any sense?).
I looked over your code a little and you're doing something a bit odd with
the predict matrix. The prediction for the CSD model should
be done from the model coefficients and not from an ODF. I'll comment more
inline to help with this issue.

Alternatively, we could implement an alternative model.predict that
takes
parameters as input and does the prediction. Is that what you were
thinking? In a way, it's similar to what we do in the sims module.

Yep that exactly what I mean, Ideally I'd like to be able to do something
like this:

model = SomeModel(gtab, ...)
sim = model.predict(params)

fit = model.fit(data)
fit.predict() == model.predict(fit.params)

MrBago · 2014-03-21T01:20:21Z

dipy/reconst/shm.py

+                self.model.cache_set("prediction_matrix",
+                                     sphere,
+                                     prediction_matrix)
+        return prediction_matrix.T


This function has a few issues. First you want to predict signal from the coefficients of the csd fit. You can do this by using simply taking dot(Pmatrix, fit.shm_coeff). In this case Pmatrix == dot(SH_basis * model.R). You can take a look at the fit method to see how we get the SH-basis and model.R is the spherical harmonic representation of the response function, meaning that coef-SIG = dot(R, coef-FOD). You can re-compute R if you want, but it's already computed for every CSD model.

Keep in mind with the response function that it needs to be axially symmetric, this is taken care of by the sh_to_rh function when R is created. We've actually done something very weird with the response function, we represent it as a tensor and that's not very good. What actually happens is that an axially symmetric SH function is fit to the tensor and that's the true response function.

arokem · 2014-06-18T18:21:47Z

Alright. This is ready for a re-review. In particular, I have implemented @MrBago's recent suggestions. Please take a look when you get a chance.

Garyfallidis · 2014-06-18T20:59:54Z

dipy/reconst/csdeconv.py

+    return pred_sig
+
+
+class ConstrainedSphericalDeconvModel(OdfModel, Cache):


Hmm... why you have moved the classes at the end of the file? Please put them up again!

This follows the general rule that things are used after they are defined.

Well... @stefanv may have a different opinion on that. Also this rule makes sense for C code but not for Python.

I personally think that when classes are defined at the top it is easier to get to the most important parts of the code faster. The CSDModel is a good example.
But because people reviewing my code have asked me in the past to put classes'
definitions at the top and I started using this style... I saw the benefits myself and
I wanted to share this experience with you too.

I've found that it's helpful to know where to look for something. If you
see something (a function, variable, etc.) being used in code, you know you
will always find it defined, allocated, or imported somewhere earlier in
the file. So, you always know which way to search.

On Wed, Jun 18, 2014 at 5:05 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

In dipy/reconst/csdeconv.py:

if np.iterable(S0):

# If it's an array, we need to give it one more dimension:

S0 = S0[..., None]

This is the key operation: convolve and multiply by S0:

pre_pred_sig = S0 * np.dot(predict_matrix, sh_coeff)

Now put everything in its right place:

pred_sig = np.zeros(pre_pred_sig.shape[:-1] + (gtab.bvals.shape[0],))

pred_sig[..., ~gtab.b0s_mask] = pre_pred_sig

pred_sig[..., gtab.b0s_mask] = S0

return pred_sig

+class ConstrainedSphericalDeconvModel(OdfModel, Cache):

Well... @stefanv https://github.com/stefanv may have a different
opinion on that. Also this rule makes sense for C code but not for Python.

I personally think that when classes are defined at the top it is easier
to get to the most important parts of the code faster. The CSDModel is a
good example.
But because people reviewing my code have asked me in the past to put
classes'
definitions at the top and I started using this style... I saw the
benefits myself and
I wanted to share this experience with you too.

—
Reply to this email directly or view it on GitHub
https://github.com/nipy/dipy/pull/335/files#r13947669.

Either way, we need to decide on a standard and stick to it everywhere. As
it is, parts of the code (e.g. reconst.dti) adhere to one convention, and
other parts adhere to another. That's certainly more confusing than any one
of these conventions by itself.

On Wed, Jun 18, 2014 at 6:48 PM, Ariel Rokem arokem@gmail.com wrote:

I've found that it's helpful to know where to look for something. If you
see something (a function, variable, etc.) being used in code, you know you
will always find it defined, allocated, or imported somewhere earlier in
the file. So, you always know which way to search.

On Wed, Jun 18, 2014 at 5:05 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

In dipy/reconst/csdeconv.py:

if np.iterable(S0):

# If it's an array, we need to give it one more dimension:

S0 = S0[..., None]

This is the key operation: convolve and multiply by S0:

pre_pred_sig = S0 * np.dot(predict_matrix, sh_coeff)

Now put everything in its right place:

pred_sig = np.zeros(pre_pred_sig.shape[:-1] + (gtab.bvals.shape[0],))

pred_sig[..., ~gtab.b0s_mask] = pre_pred_sig

pred_sig[..., gtab.b0s_mask] = S0

return pred_sig

+class ConstrainedSphericalDeconvModel(OdfModel, Cache):

Well... @stefanv https://github.com/stefanv may have a different
opinion on that. Also this rule makes sense for C code but not for Python.

I personally think that when classes are defined at the top it is easier
to get to the most important parts of the code faster. The CSDModel is a
good example.
But because people reviewing my code have asked me in the past to put
classes'
definitions at the top and I started using this style... I saw the
benefits myself and
I wanted to share this experience with you too.

—
Reply to this email directly or view it on GitHub
https://github.com/nipy/dipy/pull/335/files#r13947669.

I don't know if I shared an opinion about this in the past, but I don't see any problem with the way @arokem implemented it.

Yes - you made this comment about the dti module, when I was working on some of the metrics (that are now implemented at the top of the file). And I've followed it since. It's worked very nicely for me so far.

Ah, yes, so that makes sense to me: to define things so that any name seen must already have been defined (I somehow thought @Garyfallidis suggested that I said something different in the past).

Yes, you asked me to change the CSD class to put it on the top of the file. Don't you remember?

OK - I moved the classes back to the top of the file. We can have the
discussion about the general policy some other time.

On Wed, Jul 2, 2014 at 12:22 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

In dipy/reconst/csdeconv.py:

if np.iterable(S0):

# If it's an array, we need to give it one more dimension:

S0 = S0[..., None]

This is the key operation: convolve and multiply by S0:

pre_pred_sig = S0 * np.dot(predict_matrix, sh_coeff)

Now put everything in its right place:

pred_sig = np.zeros(pre_pred_sig.shape[:-1] + (gtab.bvals.shape[0],))

pred_sig[..., ~gtab.b0s_mask] = pre_pred_sig

pred_sig[..., gtab.b0s_mask] = S0

return pred_sig

+class ConstrainedSphericalDeconvModel(OdfModel, Cache):

Yes, you asked me to change the CSD class to put it on the top of the file.

—
Reply to this email directly or view it on GitHub
https://github.com/nipy/dipy/pull/335/files#r14477908.

arokem · 2014-06-18T23:19:27Z

Thanks for the comments @Garyfallidis. Still working on an example. I should have that up later today.

Garyfallidis · 2014-06-19T00:33:59Z

dipy/reconst/tests/test_dti.py

    assert_array_almost_equal(dmfit.predict(gtab, S0=100), S)
-
+    assert_array_almost_equal(dm.predict(dmfit.model_params, S0=100), S)


Why you selected to have the predict method at the Model class and not at the Fit class?
I would think that because you have already the params after the model.fit call then it would make sense
to call predict from the Fit. Can you elaborate a bit on this?

It's now implemented in both the Model and the Fit, following @MrBago's
suggestion. I think that's pretty nice. You can use the API you prefer,
depending on your use-case. Indeed, in the kfold cross-validation, I use
the Fit object's predict method.

On Wed, Jun 18, 2014 at 5:33 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

In dipy/reconst/tests/test_dti.py:

assert_array_almost_equal(dmfit.predict(gtab, S0=100), S)

assert_array_almost_equal(dm.predict(dmfit.model_params, S0=100), S)

Why you selected to have the predict method at the Model class and not at
the Fit class?
I would think that because you have already the params after the model.fit
call then it would make sense
to call predict from the Fit. Can you elaborate a bit on this?

—
Reply to this email directly or view it on GitHub
https://github.com/nipy/dipy/pull/335/files#r13948572.

arokem · 2014-06-20T04:02:13Z

The example is now ready for review. Thanks!

arokem · 2014-06-25T16:33:52Z

Does anyone have the chance to take a look at the example? Other than that - any more comments? @MrBago - were all your previous comments properly addressed?

stefanv · 2014-06-25T17:26:01Z

dipy/reconst/csdeconv.py

+        """
+        Predict a signal from sh coefficients for this model
+        """
+        return csd_predict(sh_coeff, self.gtab, response=self.response, S0=S0,


Is csd_predict meant to be exposed to the public? Otherwise rename _csd_predict.

I don't see any reason that this should not be publicly exposed. I have now also documented the return value.

arokem · 2014-06-27T04:57:11Z

Thanks @stefanv for that over-zealous review ;-)

I hope that I have addressed most of your comments. In particular, I would like to hear more about what you were thinking about when you wrote that comment about scikit-learn. I'd be happy to know whether they have tools that we could just somehow apply to this problem and would do this for us. I doubt that's the case. Consider the handling of gradient-tables, S0, etc.

arokem · 2014-06-27T17:03:18Z

OK - I think that I have addressed all of @stefanv's comments. At the very least, I have responded to all of them :-)

Anything more? Are we ready to merge this thing?

stefanv · 2014-06-27T17:10:12Z

@arokem Any reason in particular you think the review was over-zealous? Along with the smiley, I'm not sure whether I should tone down any specific aspect to make my comments more acceptable. Usually when I do reviews the comments are meant as recommendations, so I leave it up to the author to decide whether to respond or not (and you did, to all of them).

arokem · 2014-06-27T17:15:54Z

Sorry - didn't mean to be too facetious. I was just amused by how detailed
it was. A good thing, probably. Didn't mean to offend. Well, just poke a
bit of fun at you really. Apologize if I did. I hope that I addressed most
of the comments.

On Fri, Jun 27, 2014 at 10:10 AM, Stefan van der Walt <
notifications@github.com> wrote:

@arokem https://github.com/arokem Any reason in particular you think
the review was over-zealous? Along with the smiley, I'm not sure whether I
should tone down any specific aspect to make my comments more acceptable.
Usually when I do reviews the comments are meant as recommendations, so I
leave it up to the author to decide whether to respond or not (and you did,
to all of them).

—
Reply to this email directly or view it on GitHub
#335 (comment).

arokem · 2014-06-27T17:16:55Z

Oh sorry - still need to do round-tripping. Don't merge yet!

stefanv · 2014-06-27T17:16:55Z

Ah, no worries--I like people poking fun at me (as much as I dislike offending!). I am a total OCD programming nazi :)

arokem · 2014-06-27T17:44:35Z

OK - now it's really ready (I think)!

stefanv · 2014-06-27T17:48:10Z

Thanks, the doc update clarified things for me!

arokem · 2014-06-27T17:49:04Z

Cool. Thanks again for the careful read.

On Fri, Jun 27, 2014 at 10:48 AM, Stefan van der Walt <
notifications@github.com> wrote:

Thanks, the doc update clarified things for me!

—
Reply to this email directly or view it on GitHub
#335 (comment).

arokem · 2014-07-02T18:09:45Z

Is this ready to be merged? Any more comments?

Garyfallidis · 2014-07-02T18:31:52Z

Give me this afternoon to have a last check and if all is good will merge it.

arokem · 2014-07-02T18:40:37Z

Sweet. Thanks!

On Wed, Jul 2, 2014 at 11:31 AM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

Give me this afternoon to have a last check and if all is good will merge
it.

—
Reply to this email directly or view it on GitHub
#335 (comment).

Garyfallidis · 2014-07-02T20:39:06Z

Nice work. Now @maurozucchelli and @mpaquette can write their predict functions for the SHORE etc and test them with your framework.

Xval

arokem · 2014-07-02T20:45:25Z

Yip! GGT!

On Wed, Jul 2, 2014 at 1:39 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

Merged #335 #335.

—
Reply to this email directly or view it on GitHub
#335 (comment).

MrBago reviewed Mar 20, 2014
View reviewed changes

MrBago reviewed Mar 21, 2014
View reviewed changes

Garyfallidis reviewed Jun 18, 2014
View reviewed changes

Garyfallidis reviewed Jun 19, 2014
View reviewed changes

arokem mentioned this pull request Jun 20, 2014

CSD fit issue #340

Merged

stefanv reviewed Jun 25, 2014
View reviewed changes

arokem added 10 commits June 26, 2014 21:11

DOC: Return values of csd_predict.

da5683e

DOC: inline with mathdollar.

a86ac3c

DOC: full stop.

136bc03

RF: Better to use an explicit test.

3467cc9

DOC: Clarify.

a17ad90

DOC: Clarify.

341a79f

DOC: Explicitly state the dimensions of the SH representation.

2edf7bf

DOC: back-ticks n' stuff.

774885f

DOC: Descriptions of input parameters.

54e02d2

DOC + RF: Variable descriptions, names. A few changes in example.

628d292

RF: Changed the name of this module to be clearer.

24b9242

DOC: Clarification.

73697d8

TST: Added roundtrip testing of CSD prediction API.

a05799f

RF: Moved the class definitions back to the top of the file.

5c6bb90

Garyfallidis added a commit that referenced this pull request Jul 2, 2014

Merge pull request #335 from arokem/xval

8019449

Xval

Garyfallidis merged commit 8019449 into dipy:master Jul 2, 2014

		return pred_sig


		class ConstrainedSphericalDeconvModel(OdfModel, Cache):

		assert_array_almost_equal(dmfit.predict(gtab, S0=100), S)

		assert_array_almost_equal(dm.predict(dmfit.model_params, S0=100), S)

Xval #335

Xval #335

Conversation

arokem commented Mar 11, 2014

Garyfallidis commented Mar 17, 2014

arokem commented Mar 17, 2014

arokem commented Mar 17, 2014

MrBago commented Mar 18, 2014

arokem commented Mar 18, 2014

Choose a reason for hiding this comment

MrBago commented Mar 21, 2014

Choose a reason for hiding this comment

arokem commented Jun 18, 2014

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

This is the key operation: convolve and multiply by S0:

Now put everything in its right place:

Choose a reason for hiding this comment

This is the key operation: convolve and multiply by S0:

Now put everything in its right place:

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

This is the key operation: convolve and multiply by S0:

Now put everything in its right place:

arokem commented Jun 18, 2014

Choose a reason for hiding this comment

Choose a reason for hiding this comment

arokem commented Jun 20, 2014

arokem commented Jun 25, 2014

Choose a reason for hiding this comment

Choose a reason for hiding this comment

arokem commented Jun 27, 2014

arokem commented Jun 27, 2014

stefanv commented Jun 27, 2014

arokem commented Jun 27, 2014

arokem commented Jun 27, 2014

stefanv commented Jun 27, 2014

arokem commented Jun 27, 2014

stefanv commented Jun 27, 2014

arokem commented Jun 27, 2014

arokem commented Jul 2, 2014

Garyfallidis commented Jul 2, 2014

arokem commented Jul 2, 2014

Garyfallidis commented Jul 2, 2014

arokem commented Jul 2, 2014