Skip to content
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.

Sign up for GitHub

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

Jump to bottom

ENH: avoid calling log() on zero-valued elements in anisotropic_power #1086

Closed
wants to merge 5 commits into from
Closed

ENH: avoid calling log() on zero-valued elements in anisotropic_power #1086

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
9991a2d
DOC: some typo corrections.
jchoude Jun 22, 2016
851cb6a
ENH: only apply log on non-zero elements.
jchoude Jun 22, 2016
079d01f
PEP8
jchoude Jun 23, 2016
777dfa2
DOC: better comment.
jchoude Jun 23, 2016
6d9916a
BF: fixing 0-d indexing error on numpy 1.7.1 in tests.
jchoude Oct 26, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
20 changes: 13 additions & 7 deletions dipy/reconst/shm.py
View file
Open in desktop
Expand Up @@ -1010,7 +1010,7 @@ def anisotropic_power(sh_coeffs, norm_factor=0.00001, power=2,
----------
sh_coeffs : ndarray
A ndarray where the last dimension is the
SH coeff estimates for that voxel.
SH coefficients estimates for that voxel.
norm_factor: float, optional
The value to normalize the ap values. Default is 10^-5.
power : int, optional
Expand All @@ -1026,19 +1026,20 @@ def anisotropic_power(sh_coeffs, norm_factor=0.00001, power=2,

Notes
----------
Calculate AP image based on a IxJxKxC SH coeffecient matrix based on the
Calculate AP image based on a IxJxKxC SH coefficient matrix based on the
equation:
.. math::
AP = \sum_{l=2,4,6,...}{\frac{1}{2l+1} \sum_{m=-l}^l{|a_{l,m}|^n}}

Where the last dimension, C, is made of a flattened array of $l$x$m$
coefficients, where $l$ are the SH orders, and $m = 2l+1$,
So l=1 has 1 coeffecient, l=2 has 5, ... l=8 has 17 and so on.
A l=2 SH coeffecient matrix will then be composed of a IxJxKx6 volume.
A l=2 SH coefficient matrix will then be composed of a IxJxKx6 volume.
The power, $n$ is usually set to $n=2$.

The final AP image is then shifted by -log(normal_factor), to be strictly non-negative. Remaining values < 0 are discarded (set to 0), per default,
and this option is controlled throug the `non_negative` key word argument.
The final AP image is then shifted by -log(norm_factor), to be strictly
non-negative. Remaining values < 0 are discarded (set to 0), per default,
and this option is controlled through the `non_negative` keyword argument.

References
----------
Expand All @@ -1060,8 +1061,13 @@ def anisotropic_power(sh_coeffs, norm_factor=0.00001, power=2,
ap += ap_i
n_start = n_stop

# Shift the map to be mostly non-negative:
log_ap = np.log(ap) - np.log(norm_factor)
# Shift the map to be mostly non-negative,
# only applying the log operation to positive elements
# to avoid getting numpy warnings on log(0).
# It is impossible to get ap values smaller than 0.
# Also avoids getting voxels with -inf when non_negative=False.
log_ap = np.zeros_like(ap)
log_ap[ap > 0] = np.log(ap[ap > 0]) - np.log(norm_factor)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is more for my intellectual curiosity, does anyone knows if there is any substantial speedup of caching this boolean indexing (i.e. ap > 0) or does Python will optimize it anyway?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't know whether Python has some optimization in place, but there seems to be only a small advantage to caching:
https://gist.github.com/arokem/fba319fedf2b541240504f56f97af45e

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had made quick tests to see if it was worth it, and I also didn't see a significant improvement. However, I'm no expert on the internals of Numpy or python, so there might be some magic caching behind the scenes, or an optimized accessor.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, good to know, thanks.


# Deal with residual negative values:
if non_negative:
Expand Down
2 changes: 1 addition & 1 deletion dipy/reconst/tests/test_shm.py
View file
Open in desktop
Expand Up @@ -437,7 +437,7 @@ def test_anisotropic_power():
assert_array_almost_equal(apvals, answers)
# Test that this works for single voxel arrays as well:
assert_array_almost_equal(
anisotropic_power(coeffs[1], norm_factor=norm_factor),
anisotropic_power(np.asarray(coeffs[1]), norm_factor=norm_factor),
answers[1])


Expand Down