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

Faster spherical harmonics #413

Merged
merged 6 commits into from Sep 16, 2014

Conversation

Projects
None yet
4 participants
@Garyfallidis
Member

Garyfallidis commented Sep 5, 2014

Thanx @ChantalTax, @pv and @stefanv for helping with this. Here is a faster implementation of spherical harmonics. This should make the calculation of real spherical harmonics and complex spherical harmonics faster everywhere in dipy.

Let me know what you think! This is now 100X faster. :D

def spherical_harmonics(m, n, theta, phi):
r""" Compute spherical harmonics
This may take scalar or array arguments. The inputs will be broadcasted

This comment has been minimized.

@arokem

arokem Sep 5, 2014

Member

This all looks awesome to me. My only comment is an English thing: I believe this should be simply broadcast

http://www.oed.com/view/Entry/23508?redirectedFrom=broadcasted#eid13325489

Also - what exactly do you mean by "inputs will be broadcast against each other"? theta and phi will be made to conform if they can be broadcast to each other's shape?

This comment has been minimized.

@Garyfallidis

Garyfallidis Sep 5, 2014

Member

I used what it was in the documentation of scipy.special.sph_harm
http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.special.sph_harm.html

I think it is fine to say broadcasted. Even if it is not in some dictionaries. It is commonly used.

@@ -171,6 +206,9 @@ def real_sph_harm(m, n, theta, phi):
The azimuthal (longitudinal) coordinate.
phi : float [0, pi]
The polar (colatitudinal) coordinate.
fast : bool
If False use the scipy.special.sph_harm else use spherical_harmonics
from dipy.reconst.shm.

This comment has been minimized.

@MrBago

MrBago Sep 5, 2014

Contributor

Why are we doing this? Why give the user the option? If we want to give the user the option this doesn't seem like an appropriate way to do it.

This comment has been minimized.

@arokem

arokem Sep 6, 2014

Member

But sure is useful for regression testing :-)

On Fri, Sep 5, 2014 at 2:53 PM, MrBago notifications@github.com wrote:

In dipy/reconst/shm.py:

@@ -171,6 +206,9 @@ def real_sph_harm(m, n, theta, phi):
The azimuthal (longitudinal) coordinate.
phi : float [0, pi]
The polar (colatitudinal) coordinate.

  • fast : bool
  •    If False use the scipy.special.sph_harm else use spherical_harmonics
    
  •    from dipy.reconst.shm.
    

Why are we doing this? Why give the user the option? If we want to give
the user the option this doesn't seem like an appropriate way to do it.


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

This comment has been minimized.

@pv

pv Sep 6, 2014

A faster sph_harm was added to Scipy in May.

This comment has been minimized.

@arokem

arokem Sep 6, 2014

Member

Another option for handling this is to separate out the new implementation
of the sph_harm, and depending on the version of scipy that is installed,
either import from scipy or from our own. That might make ultimately
deprecating this, as the scipy implementation percolates through, easier
(?).

On Sat, Sep 6, 2014 at 9:41 AM, Pauli Virtanen notifications@github.com
wrote:

In dipy/reconst/shm.py:

@@ -171,6 +206,9 @@ def real_sph_harm(m, n, theta, phi):
The azimuthal (longitudinal) coordinate.
phi : float [0, pi]
The polar (colatitudinal) coordinate.

  • fast : bool
  •    If False use the scipy.special.sph_harm else use spherical_harmonics
    
  •    from dipy.reconst.shm.
    

A faster sph_harm was added to Scipy in May.


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

-1.06955885, -2.83826831, 1.81932195, 2.81296654])
rsh = real_sph_harm(m, n, theta[:, None], phi[:, None], fast=True)
rsh2 = real_sph_harm(m, n, theta[:, None], phi[:, None], fast=False)

This comment has been minimized.

@MrBago

MrBago Sep 5, 2014

Contributor

This test looks fine, but why not add a simple test to also test that scipy.sph_harm(...) == dipy.sph_harm(...).

@MrBago

This comment has been minimized.

Contributor

MrBago commented Sep 5, 2014

Any thoughts about pushing this upstream to scipy? Is there any reason that the larger scipy community would not want to benefit from this awesomeness? Or is this already in the latest scipy and I don't know about it?

@arokem

This comment has been minimized.

Member

arokem commented Sep 6, 2014

BTW - any idea why Travis is seg-faulting on this one? Seems like just a
glitch - I don't get that on my machine.

On Fri, Sep 5, 2014 at 5:01 PM, Ariel Rokem arokem@gmail.com wrote:

Seems like that's already happening. See @stefanv's comments on gh-404.

On Fri, Sep 5, 2014 at 3:01 PM, MrBago notifications@github.com wrote:

Any thoughts about pushing this upstream to scipy? Is there any reason
that the larger scipy community would not want to benefit from this
awesomeness?


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

@arokem

This comment has been minimized.

Member

arokem commented Sep 6, 2014

Seems like that's already happening. See @stefanv's comments on gh-404.

On Fri, Sep 5, 2014 at 3:01 PM, MrBago notifications@github.com wrote:

Any thoughts about pushing this upstream to scipy? Is there any reason
that the larger scipy community would not want to benefit from this
awesomeness?


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

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Sep 12, 2014

Hey guys,

With the last commits I now address all your comments. I test sph_harm and in spherical_harmonics I check if scipy version is higher that 0.15.0.

Enjoy the speed! Merge if you think is good so we can test the recalibration with real data!
:)

@@ -37,6 +37,7 @@
from dipy.core.onetime import auto_attr
from dipy.reconst.cache import Cache
import scipy.version as SCIPY_VERSION

This comment has been minimized.

@MrBago

MrBago Sep 12, 2014

Contributor

I think this should be scipy.version.version.

"""
Compute real spherical harmonics, where the real harmonic $Y^m_n$ is
defined to be:
if str(SCIPY_VERSION) >= '0.15.0':

This comment has been minimized.

@MrBago

MrBago Sep 12, 2014

Contributor

I think "0.9.0" > "0.15.0", try using distutils.version.StrictVersion.

Compute real spherical harmonics, where the real harmonic $Y^m_n$ is
defined to be:
if str(SCIPY_VERSION) >= '0.15.0':
return sph_harm(m, n, theta, phi)

This comment has been minimized.

@MrBago

MrBago Sep 12, 2014

Contributor

Why not put the check at the module level so it's not done every time the function is called, ie:

import scipy import foo
def bar():
    pass
if scipy_version > "0.15":
    foo = bar
@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Sep 12, 2014

Ok thx @MrBago! Ready here!

@Garyfallidis

This comment has been minimized.

Member

Garyfallidis commented Sep 16, 2014

What are we waiting for here? Can we merge this and then finalize Chantal's PR too. @MrBago, @arokem ?

arokem added a commit that referenced this pull request Sep 16, 2014

@arokem arokem merged commit 924519a into nipy:master Sep 16, 2014

1 check passed

continuous-integration/travis-ci The Travis CI build passed
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment