New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Sparse Fascicle Model #460

Merged
merged 71 commits into from Dec 18, 2014

Conversation

Projects
None yet
6 participants
@arokem
Member

arokem commented Nov 11, 2014

This is a start on an implementation of the model we describe here:

http://arxiv.org/abs/1411.0721

The main current issue with this is that this will have a local dependency on scikit-learn, because it uses ElasticNet for the fit. How do people feel about that?

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 11, 2014

Option skimage dependency is no problem IMO; is that possible?

@arokem

This comment has been minimized.

Member

arokem commented Nov 11, 2014

Not following: you think we could make skimage an optional dependency, instead of sklearn? Or did you actually mean sklearn? Certainly: optional, either way.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 11, 2014

Sorry - yes - sklearn - I think it's not optional at the moment?

You could use optpkg for this, example:

./dipy/core/profile.py:
    7 : from ..utils.optpkg import optional_package
@arokem

This comment has been minimized.

Member

arokem commented Nov 11, 2014

Thanks - very helpful. How do you feel about falling back on
scipy.optimize.nnls when sklearn is not available? It would still work, but
probably not as well. Or is it preferable to just raise an error, and ask
the user to go install sklearn?

On Tue, Nov 11, 2014 at 2:04 PM, Matthew Brett notifications@github.com
wrote:

Sorry - yes - sklearn - I think it's not optional at the moment?

You could use optpkg for this, example:

./dipy/core/profile.py:
7 : from ..utils.optpkg import optional_package


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

@arokem arokem changed the title from WIP: Starting to implement the SFM. to Sparse Fascicle Model Nov 15, 2014

@arokem

This comment has been minimized.

Member

arokem commented Nov 16, 2014

I believe this is now ready for a review! @stefanv - you know this stuff pretty well. Would you mind taking a look? Thanks!

@arokem arokem force-pushed the arokem:sparse-fascicle branch from f2ef9fe to f77a466 Nov 17, 2014

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 18, 2014

I'm not sure; do you think the sklearn and scipy versions would give very different results? If so, that might be a bit strange.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 18, 2014

I mean, that would puzzling for people comparing results across machines and installs.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 18, 2014

Why wouldn't they, non negative least square and elastic net are two different solver, the first one using a positivity constraint on the solution and the second one a L1 and L2 regularization.

https://en.wikipedia.org/wiki/Elastic_net_regularization
https://en.wikipedia.org/wiki/Non-negative_least_squares

I'd actually be surprised if they gave similar results. To ensure consistency it might be better to choose one or the other instead of a fallback.

@arokem

This comment has been minimized.

Member

arokem commented Nov 18, 2014

That is true. Notice that ElasticNet is also used with a non-negativity
constraint. But, the solutions can be quite different. That said, from my
own experiments, NNLS and ElasticNet can have rather similar
cross-validation quality. Here's one such experiment:

http://nbviewer.ipython.org/gist/arokem/4253775a276d63852c7d

The bottom graph shows a comparison of cross-validation error (relative
RMSE, smaller numbers are better) for NNLS and ElasticNet. Almost as good,
even with quite different FODs...

I am thinking that instead of making this an automatic fall-back, we can
make this an input where:

solver='ElasticNet'

Is the default, and

solver='NNLS'

Uses opt.nnls instead.

What do you think?

On Mon, Nov 17, 2014 at 5:26 PM, Samuel St-Jean notifications@github.com
wrote:

Why would they, non negative least square and elastic net are two
different solver, the first one using a positivity constraint on the
solution and the second one a L1 and L2 regularization.

https://en.wikipedia.org/wiki/Elastic_net_regularization
https://en.wikipedia.org/wiki/Non-negative_least_squares

I'd actually be surprised if they gave similar results. To ensure
consistency it might be better to choose one or the other instead of a
fallback.


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

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 18, 2014

Having the default be explicit sounds like a good solution to the problem.

@samuelstjean

This comment has been minimized.

Contributor

samuelstjean commented Nov 18, 2014

Wouldn't having the fallback work less good than the original option instead cripple the function? Instead of having something somewhat working, might be better to force the use of sklearn instead of having suboptimal results.

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 18, 2014

I think the idea is that the default is the sklearn elasticnet optimization, and the routine will error (by default) if sklearn is not installed, rather than falling back to NNLS. But, you can pass in an explicit flag to use the NNLS optimization.

@arokem

This comment has been minimized.

Member

arokem commented Nov 19, 2014

I've tried to go with a more general approach. We are working on some other algorithms to perform this fitting, so providing an interface to this is not just a theoretical exercise.

I can't seem to get the implementation of the testing for membership in a super-class quite right (here: https://github.com/arokem/dipy/blob/beb10c4130baa089098a41f7bfa87e9ab9f62d24/dipy/reconst/sfm.py#L147). I am not sure how this is done, but I'm sure that there is some standard way. Unfortunately, I don't have time right this moment to look into that. I hope to get back to that some time later this week

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 19, 2014

isinstance not what you want there?

@arokem

This comment has been minimized.

Member

arokem commented Nov 20, 2014

It seems not (see most recent commit).

When dropping into a debugger, I get:


(Pdb) type(solver)

Am I doing the inheritance wrong?

On Wed, Nov 19, 2014 at 2:08 PM, Matthew Brett notifications@github.com
wrote:

isinstance not what you want there?


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

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Nov 20, 2014

Wow, that is confusing. Is it easy for me to replicate what you did?

@arokem

This comment has been minimized.

Member

arokem commented Nov 20, 2014

OK - I figured out what that was all about. I was passing in SillySolver,
instead of SillySolver(). The class, instead of an instance ...

On Wed, Nov 19, 2014 at 5:10 PM, Matthew Brett notifications@github.com
wrote:

Wow, that is confusing. Is it easy for me to replicate what you did?


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

@arokem

This comment has been minimized.

Member

arokem commented Nov 20, 2014

What do you think about the general approach?

On Wed, Nov 19, 2014 at 5:27 PM, Ariel Rokem arokem@gmail.com wrote:

OK - I figured out what that was all about. I was passing in
SillySolver, instead of SillySolver(). The class, instead of an
instance ...

On Wed, Nov 19, 2014 at 5:10 PM, Matthew Brett notifications@github.com
wrote:

Wow, that is confusing. Is it easy for me to replicate what you did?


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

@@ -237,3 +240,65 @@ def __del__(self):
for fname in self.tmp_files:
os.remove(fname)
class SKLearnLinearSolver(with_metaclass(abc.ABCMeta, object)):

This comment has been minimized.

@stefanv

stefanv Nov 21, 2014

Contributor

It's probably not necessary to add the , object part here.

Provide a sklearn-like uniform interface to algorithms that solve problems
of the form:
.. math::

This comment has been minimized.

@stefanv

stefanv Nov 21, 2014

Contributor

I think we have the math-dollar extension active, so can simply be $y = Ax$ and $x$.

Sub-classes of SKLearnLinearSolver should provide a 'fit' method that have
the following signature: `SKLearnLinearSolver.fit(X, y)`, which would set
an attribute `SKLearnLinearSolver.coef_`, with the shape (X.shape[1],),
such that an estimate of y can be calculated as: SKLearnLinearSolver

This comment has been minimized.

@stefanv

stefanv Nov 21, 2014

Contributor

can be calculated as... ?

"""Implement for all derived classes """
def predict(self, X):

This comment has been minimized.

@stefanv

stefanv Nov 21, 2014

Contributor

Why not also make this an abstract method?

C : array, shape = (n_samples,)
Predicted values.
"""
X = np.asarray(X)

This comment has been minimized.

@stefanv

stefanv Nov 21, 2014

Contributor

Why the custom implementation, instead of using the scikit-learn solver's predict method?

This comment has been minimized.

@arokem

arokem Nov 24, 2014

Member

It requires scikit-learn, and we're trying to deal with situations where that's not installed. Though, why wouldn't it?

The eigenvalues of a tensor which will serve as a kernel function
mode : str
'sig' : for a signal design matrix. 'odf' for an odf convolution matrix

This comment has been minimized.

@stefanv

stefanv Nov 21, 2014

Contributor

I'd expand 'sig' to 'signature' and expand this explanation a bit.

This comment has been minimized.

@arokem

arokem Nov 24, 2014

Member

Expanded this to the full signal. Is that what you meant?

Sets the columns of the matrix
response : list of 3 elements
The eigenvalues of a tensor which will serve as a kernel function

This comment has been minimized.

@stefanv

stefanv Nov 21, 2014

Contributor

An example below may help the user to pick reasonable values.

Returns
-------
mat : ndarray
A matrix either for deconvolution with the signal, or for reconvolution

This comment has been minimized.

@stefanv

stefanv Nov 21, 2014

Contributor

A little bit sparse--perhaps expand this sentence as well?

This comment has been minimized.

@stefanv

stefanv Nov 25, 2014

Contributor

How about this one?

# predicted by a "canonical", symmetrical tensor rotated towards this
# vertex of the sphere:
canonical_tensor = np.array([[response[0], 0, 0],
[0, response[1], 0],

This comment has been minimized.

@stefanv

stefanv Nov 21, 2014

Contributor

indentation

This comment has been minimized.

@stefanv

stefanv Nov 21, 2014

Contributor

np.diag(response)

mat = np.empty((gtab.x.shape[0], sphere.vertices.shape[0]))
# Calculate column-wise:
for ii, this_dir in enumerate(sphere.vertices):

This comment has been minimized.

@stefanv

stefanv Nov 21, 2014

Contributor

Nicely done

@arokem

This comment has been minimized.

Member

arokem commented Dec 16, 2014

I've refactored the examples quite thoroughly, using the new tracking API, etc. Please take a look and let me know whether it looks OK.

"""
Next, we will create a visualization of these streamlines, relative to this
subjects T1-weighted anatomy:

This comment has been minimized.

@stefanv

stefanv Dec 16, 2014

Contributor

subject's

@stefanv

This comment has been minimized.

Contributor

stefanv commented Dec 16, 2014

Beautiful example.

@arokem

This comment has been minimized.

Member

arokem commented Dec 16, 2014

Thanks for the feedback. Typo fixed!

On Tue, Dec 16, 2014 at 3:44 AM, Stefan van der Walt <
notifications@github.com> wrote:

Beautiful example.


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

sphere,
relative_peak_threshold=.5,
min_separation_angle=25,
return_sh=False)

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 16, 2014

Member

You forgot parallel=True

sufficiently high.
"""
from dipy.tracking.local import ThresholdTissueClassifier

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 16, 2014

Member

Hi, here you changed your previous tutorial which was full brain tracking and copied some parts from the new tutorial from the new tracking API. Now we have two very similar tutorials on tracking with only the model being different.
What if you use full brain tracking here and not seed from the CC? This will have an added value because you show something that it is not available in the other tutorials.

.. _sfm-track:
==================================================
Tracking with the Sparse Fascicle Model and EuDX

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 16, 2014

Member

Here you need to change the title. You don't use EuDX anymore.

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Dec 16, 2014

Hi @arokem, indeed great work as I said before. I gave a suggestion for the tracking example. If you think you want to work on that I would suggest to do it on a separate PR and remove the tracking example from this PR. Apart from this. This is ready to be merged.

Here is also a feature request. I would like to see in the future if possible the actual tensors coming out of the SFM. I think this would be very useful as you could start with those solutions and then refine them to do a complete MultiTensor fitting which would need much more time otherwise. Interesting?

@arokem

This comment has been minimized.

Member

arokem commented Dec 16, 2014

Cool. I just made some changes to the tracking example, as suggested.
Running it right now, and will push into that same PR, when I am done
running it.

On Tue, Dec 16, 2014 at 12:18 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

Hi @arokem https://github.com/arokem, indeed great work as I said
before. I gave a suggestion for the tracking example. If you think you want
to work on that I would suggest to do it on a separate PR and remove the
tracking example from this PR. Apart from this. This is ready to be merged.

Here is also a feature request. I would like to see in the future if
possible the actual tensors coming out of the SFM. I think this would be
very useful as you could start with those solutions and then refine them to
do a complete MultiTensor fitting which would need much more time
otherwise. Interesting?


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

@arokem

This comment has been minimized.

Member

arokem commented Dec 16, 2014

Concerning the feature request: I do think it's interesting, and we've
already started some work in that direction, so you can hope to see that in
some future version of this.

On Tue, Dec 16, 2014 at 12:23 PM, Ariel Rokem arokem@gmail.com wrote:

Cool. I just made some changes to the tracking example, as suggested.
Running it right now, and will push into that same PR, when I am done
running it.

On Tue, Dec 16, 2014 at 12:18 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

Hi @arokem https://github.com/arokem, indeed great work as I said
before. I gave a suggestion for the tracking example. If you think you want
to work on that I would suggest to do it on a separate PR and remove the
tracking example from this PR. Apart from this. This is ready to be merged.

Here is also a feature request. I would like to see in the future if
possible the actual tensors coming out of the SFM. I think this would be
very useful as you could start with those solutions and then refine them to
do a complete MultiTensor fitting which would need much more time
otherwise. Interesting?


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

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Dec 16, 2014

Super!

@yarikoptic

This comment has been minimized.

Member

yarikoptic commented Dec 16, 2014

sorry -- I have missed that there was a question to me -- sure, drop 2.6 if it helps! even current debian stable is 2.7 ;)

@arokem

This comment has been minimized.

Member

arokem commented Dec 16, 2014

Thanks @yarikoptic for that info. For now, it seems that we keep 2.6, see:
#471
But we can drop it eventually

arokem added some commits Dec 16, 2014

from dipy.tracking.local import LocalTracking
streamlines = LocalTracking(pnm, classifier, seeds, affine, step_size=.5)
streamlines = list(streamlines)

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 17, 2014

Member

If someone does here what you suggest and remove that line it will crash his memory because later in fvtk.streamtube the program will try to create streamtubes for too many streamlines. What I would do is to use all the seeds here but later in the visualization I would use select_random_set_of_streamlines from dipy.tracking.streamline and visualize only a random set of all the streamlines generated. If you agree with the suggestion and change this I will merge this PR. Thx again.

This comment has been minimized.

@arokem

arokem Dec 17, 2014

Member

Will do ASAP.

On Tue, Dec 16, 2014 at 4:48 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

In doc/examples/sfm_tracking.py
#460 (diff):

+For the sake of brevity, we will take only the first 1000 seeds, generating
+only 1000 streamlines. Remove this line to track from many more points in all of
+the white matter
+"""
+
+seeds = seeds[:1000]
+
+"""
+We now have the necessary components to construct a tracking pipeline and
+execute the tracking
+"""
+
+from dipy.tracking.local import LocalTracking
+streamlines = LocalTracking(pnm, classifier, seeds, affine, step_size=.5)
+
+streamlines = list(streamlines)

If someone does here what you suggest and remove that line it will crash
his memory because later in fvtk.streamtube the program will try to create
streamtubes for too many streamlines. What I would do is to use all the
seeds here but later in the visualization I would use
select_random_set_of_streamlines from dipy.tracking.streamline and
visualize only a random set of all the streamlines generated. If you agree
with the suggestion and change this I will merge this PR. Thx again.


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

This comment has been minimized.

@arokem

arokem Dec 17, 2014

Member

OK - all set - I think. Do it!

On Tue, Dec 16, 2014 at 4:49 PM, Ariel Rokem arokem@gmail.com wrote:

Will do ASAP.

On Tue, Dec 16, 2014 at 4:48 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

In doc/examples/sfm_tracking.py
#460 (diff):

+For the sake of brevity, we will take only the first 1000 seeds, generating
+only 1000 streamlines. Remove this line to track from many more points in all of
+the white matter
+"""
+
+seeds = seeds[:1000]
+
+"""
+We now have the necessary components to construct a tracking pipeline and
+execute the tracking
+"""
+
+from dipy.tracking.local import LocalTracking
+streamlines = LocalTracking(pnm, classifier, seeds, affine, step_size=.5)
+
+streamlines = list(streamlines)

If someone does here what you suggest and remove that line it will crash
his memory because later in fvtk.streamtube the program will try to create
streamtubes for too many streamlines. What I would do is to use all the
seeds here but later in the visualization I would use
select_random_set_of_streamlines from dipy.tracking.streamline and
visualize only a random set of all the streamlines generated. If you agree
with the suggestion and change this I will merge this PR. Thx again.


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

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Dec 18, 2014

Sorry for the delay I had to sleep down the jet lag.

Many congratulations Ariel. Another awesome new toy for Dipy :)

Garyfallidis added a commit that referenced this pull request Dec 18, 2014

@Garyfallidis Garyfallidis merged commit 5f7512f into nipy:master Dec 18, 2014

1 check passed

continuous-integration/travis-ci The Travis CI build passed
Details
@stefanv

This comment has been minimized.

Contributor

stefanv commented Dec 18, 2014

This is awesome indeed! Thank you @arokem

@arokem

This comment has been minimized.

Member

arokem commented Dec 19, 2014

Thanks! Just updated the website with these new examples:
http://nipy.org/dipy/examples_built/sfm_tracking.html#example-sfm-tracking

On Thu, Dec 18, 2014 at 3:38 PM, Stefan van der Walt <
notifications@github.com> wrote:

This is awesome indeed! Thank you @arokem https://github.com/arokem


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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment