BF: Keep up with changes in scipy 0.16

[arokem]

Recent changes in scipy use a different heuristic to find the epsilon
parameter used for RBF interpolation, causing for us issue #688. This
should fix it, by using the old heuristic.

[arokem]

Note - the first commit (arokem@7c4e670) are the substantial changes. The rest is just PEP8 fixes that I did while I was in this file already.

[arokem]

Oh dear. Hey @stefanv - any ideas here? Looks pretty bad.

[arokem]

Just to be clear on this one - these tests would have failed with previous versions of scipy as well (this PR basically just patches this up to the previous behavior, but adds more tests). It seems that RBF is just inaccurate in many use-cases, which we weren't thinking of before.

[Garyfallidis]

Oy! This is strange. Does that mean that all the functionality with interpolating on spheres is now uncertain?

[arokem]

Well - it was always uncertain. It's just that we now know for certain just
how uncertain it was...

On Sat, Sep 12, 2015 at 4:39 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

Oy! This is strange. Does that mean that all the functionality with
interpolating on sphere is now uncertain?

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

[Garyfallidis]

Okay so what is your suggestion. How to move forward with this?

[arokem]

A few options (in order of preference, I think):

• Spend some time understanding RBF interpolation, and fix it up (might
  need to fix the scipy implementation, see links to discussions there - IIUC
  that epsilon parameter is pretty much made up there as well...).
• Put a big fat warning on RBF for users to check the results of
  interpolation when they use this, and relax test requirements.
• Drop RBF interpolation altogether.

I can spend some time on option 1, but not before this Thursday. Even then,
progress might be slow. What do you think?

On Sat, Sep 12, 2015 at 4:44 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

Okay so what is your suggestion. How to move forward with this?

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

[Garyfallidis]

Hi @arokem,
here is some feedback from my investigations.

The problem at least at my scipy 0.14 seems to be only with norm="euclidean_norm"
(norm="angle" works fine) which can be resolved using a small smoothing (here using 0.1).

Check my code snippet here

def test_interp_rbf():
    def data_func(s, a, b):
        return a * np.cos(s.theta) + b * np.sin(s.phi)

    from dipy.core.sphere import Sphere, interp_rbf
    from dipy.core.subdivide_octahedron import create_unit_hemisphere
    import numpy as np

    s0 = create_unit_hemisphere(5)
    s1 = create_unit_hemisphere(6)

    from dipy.viz import fvtk
    s0 = Sphere(xyz = np.vstack((s0.vertices, -s0.vertices)))
    s1 = Sphere(xyz = np.vstack((s1.vertices, -s1.vertices)))


    for a, b in zip([1, 2, 0.5], [1, 0.5, 2]):
        ren = fvtk.ren()
        print(a, b)
        data = data_func(s0, a, b)
        expected = data_func(s1, a, b)
        interp_data_en = interp_rbf(data, s0, s1, smooth=.1, norm="euclidean_norm")
        interp_data_a = interp_rbf(data, s0, s1, smooth=.1, norm="angle")

        sf = np.zeros((3, s1.vertices.shape[0]))
        sf[0] = expected
        sf[1] = interp_data_en
        sf[2] = interp_data_a
        fvtk.add(ren, fvtk.sphere_funcs(sf, sphere=s1, norm=False, radial_scale=False))
        nt.assert_(np.mean(np.abs(interp_data_en - expected)) < 0.1)
        nt.assert_(np.mean(np.abs(interp_data_a - expected)) < 0.1)

        fvtk.show(ren)


if __name__ == "__main__":
    # nt.run_module_suite()
    test_interp_rbf()

Now if you change to the default smooth (zero) you get some weird result looking

Some other important points to improving the tests. I think those tests should use higher sampled spheres. We were using create_unit_hemisphere(2) and (3) which provide only 9 and 33 points respectively (too discrete). Also you will need to make sure that using hemispheres is a good idea for this type of interpolation. I am afraid that the points in the edges might create problems. So, I would use full spheres to be on the safe side.

Finally, I noticed some another minor issue. In lines 550 and 552 we use is for comparing strings which is generally not recommended. I would change that to ==.

So, in summary it seems that the problem may not be so serious and it can be resolved with a tiny bit of smoothing and less discrete spheres. So, if you can update the test I think we can move forward. At the same time it would be nice to have in our mind any new developments from the scipy devs. And give them feedback as we learn more from the behaviour of these functions.

How does that sound?

[arokem]

So what do you think? Should we enforce smoothing per default? Might be a
safer way to go(?) And drop the Euclidean distance? In a way, it's not
surprising that doesn't work as well as angle - Euclidean is just an
approximation of angle in this case.

I've made the following changes:

• Smoothing is set to 0.1 per default.
• norm is set to angle per default. If a user sets it to
  euclidean_norm a warning is issued. The euclidean will be completely
  removed as an option in a subsequent release cycle. I've added a test to
  see that this warning is raised.
• I've set the test to use spheres, instead of hemispheres (good call), and
  these have more points than before.

On Sun, Sep 13, 2015 at 4:20 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

Hi @arokem https://github.com/arokem,
here is some feedback from my investigations.

The problem at least at my scipy 0.14 seems to be only with
norm="euclidean_norm"
(norm="angle" works fine) which can be resolved using a small smoothing
(here using 0.1).

Check my code snippet here

def test_interp_rbf():
def data_func(s, a, b):
return a * np.cos(s.theta) + b * np.sin(s.phi)

from dipy.core.sphere import Sphere, interp_rbf
from dipy.core.subdivide_octahedron import create_unit_hemisphere
import numpy as np

s0 = create_unit_hemisphere(5)
s1 = create_unit_hemisphere(6)

from dipy.viz import fvtk
s0 = Sphere(xyz = np.vstack((s0.vertices, -s0.vertices)))
s1 = Sphere(xyz = np.vstack((s1.vertices, -s1.vertices)))


for a, b in zip([1, 2, 0.5], [1, 0.5, 2]):
    ren = fvtk.ren()
    print(a, b)
    data = data_func(s0, a, b)
    expected = data_func(s1, a, b)
    interp_data_en = interp_rbf(data, s0, s1, smooth=.1, norm="euclidean_norm")
    interp_data_a = interp_rbf(data, s0, s1, smooth=.1, norm="angle")

    sf = np.zeros((3, s1.vertices.shape[0]))
    sf[0] = expected
    sf[1] = interp_data_en
    sf[2] = interp_data_a
    fvtk.add(ren, fvtk.sphere_funcs(sf, sphere=s1, norm=False, radial_scale=False))
    nt.assert_(np.mean(np.abs(interp_data_en - expected)) < 0.1)
    nt.assert_(np.mean(np.abs(interp_data_a - expected)) < 0.1)

    fvtk.show(ren)

if name == "main":
# nt.run_module_suite()
test_interp_rbf()

[image: image]
https://cloud.githubusercontent.com/assets/134276/9839458/23f7d4f8-5a49-11e5-8b7f-72da9618fd78.png

Now if you change to the default smooth (zero) you get some weird result
looking

[image: image]
https://cloud.githubusercontent.com/assets/134276/9839503/84ae56c2-5a4a-11e5-8615-eefafd8dd7ed.png

Some other important points to improving the tests. I think those tests
should use higher sampled spheres. We were using create_unit_hemisphere(2)
and (3) which provide only 9 and 33 points respectively. Also you will need
to make sure that using hemispheres is a good idea for this type of
interpolation. I am afraid that the points in the edges might create
problems. So, I would use full spheres to be on the safe side.

Finally, I noticed some another minor issue. In lines 550 and 552 we use
is for comparing strings which is generally not recommended. I would
change that to ==.

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

[arokem]

Thanks for all the feedback, by they way!

On Sun, Sep 13, 2015 at 7:20 PM, Ariel Rokem arokem@gmail.com wrote:

So what do you think? Should we enforce smoothing per default? Might be a
safer way to go(?) And drop the Euclidean distance? In a way, it's not
surprising that doesn't work as well as angle - Euclidean is just an
approximation of angle in this case.

I've made the following changes:

• Smoothing is set to 0.1 per default.
• norm is set to angle per default. If a user sets it to
  euclidean_norm a warning is issued. The euclidean will be completely
  removed as an option in a subsequent release cycle. I've added a test to
  see that this warning is raised.
• I've set the test to use spheres, instead of hemispheres (good call),
  and these have more points than before.

On Sun, Sep 13, 2015 at 4:20 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

Hi @arokem https://github.com/arokem,
here is some feedback from my investigations.

The problem at least at my scipy 0.14 seems to be only with
norm="euclidean_norm"
(norm="angle" works fine) which can be resolved using a small smoothing
(here using 0.1).

Check my code snippet here

def test_interp_rbf():
def data_func(s, a, b):
return a * np.cos(s.theta) + b * np.sin(s.phi)

from dipy.core.sphere import Sphere, interp_rbf
from dipy.core.subdivide_octahedron import create_unit_hemisphere
import numpy as np

s0 = create_unit_hemisphere(5)
s1 = create_unit_hemisphere(6)

from dipy.viz import fvtk
s0 = Sphere(xyz = np.vstack((s0.vertices, -s0.vertices)))
s1 = Sphere(xyz = np.vstack((s1.vertices, -s1.vertices)))


for a, b in zip([1, 2, 0.5], [1, 0.5, 2]):
    ren = fvtk.ren()
    print(a, b)
    data = data_func(s0, a, b)
    expected = data_func(s1, a, b)
    interp_data_en = interp_rbf(data, s0, s1, smooth=.1, norm="euclidean_norm")
    interp_data_a = interp_rbf(data, s0, s1, smooth=.1, norm="angle")

    sf = np.zeros((3, s1.vertices.shape[0]))
    sf[0] = expected
    sf[1] = interp_data_en
    sf[2] = interp_data_a
    fvtk.add(ren, fvtk.sphere_funcs(sf, sphere=s1, norm=False, radial_scale=False))
    nt.assert_(np.mean(np.abs(interp_data_en - expected)) < 0.1)
    nt.assert_(np.mean(np.abs(interp_data_a - expected)) < 0.1)

    fvtk.show(ren)

if name == "main":
# nt.run_module_suite()
test_interp_rbf()

[image: image]
https://cloud.githubusercontent.com/assets/134276/9839458/23f7d4f8-5a49-11e5-8b7f-72da9618fd78.png

Now if you change to the default smooth (zero) you get some weird result
looking

[image: image]
https://cloud.githubusercontent.com/assets/134276/9839503/84ae56c2-5a4a-11e5-8615-eefafd8dd7ed.png

Some other important points to improving the tests. I think those tests
should use higher sampled spheres. We were using create_unit_hemisphere(2)
and (3) which provide only 9 and 33 points respectively. Also you will need
to make sure that using hemispheres is a good idea for this type of
interpolation. I am afraid that the points in the edges might create
problems. So, I would use full spheres to be on the safe side.

Finally, I noticed some another minor issue. In lines 550 and 552 we use
is for comparing strings which is generally not recommended. I would
change that to ==.

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

[Garyfallidis]

Great! All good here.
I think this is good as it is and it doesn't affect any other parts of the code.

But put in your TODO list (after Thursday) to update the api_changes.rst file by saying that the default parameter has changed and some norms will be deprecated. You should keep an issue open for scipy 0.16 if necessary.

I will go ahead and merge this as it is now. That will make it easier to know which of the other PRs are failing for a good reason and not because of RBF.

Oh and best of luck with your deadline. Thank you for looking into this.

[stefanv]

I guess it's not too surprising that the Euclidean norm performed badly. This looks like a reasonable solution!