ENH: special: Add vectorized spherical Bessel functions.

[tpudlik]

This PR provides four new ufuncs,

• spherical_jn
• spherical_yn
• spherical_in
• spherical_kn

They are vectorized alternatives to the specfun wrappers sph_jn, sph_yn, sph_in, and sph_kn, providing a nicer interface, slightly higher performance and (thanks to the use of different algorithms) better stability for large real inputs. Derivatives of spherical Bessel functions, which are also computed by the specfun wrappers, are not yet implemented.

The performance improvement is not dramatic, at a little over an order of magnitude:

@np.vectorize
def old_way_kn(n, z):
    return sph_kn(n, z)[0][-1]

x = 10**np.linspace(-2, 3, 5000)
n = 15

>>> %timeit old_way_kn(n, x)
10 loops, best of 3: 72.4 ms per loop

>>> %timeit spherical_kn(n, x)
100 loops, best of 3: 2.57 ms per loop

Better stability, especially in spherical_jn and spherical_in for large (> 10**3) real argument values, is obtained by replacing the recursive specfun algorithms with calls to AMOS/CEPHES cylindrical Bessel function routines, using these identities.

Fixes #1641 and #2165, and (except for the still missing derivatives) implements #3113.

[tpudlik]

Derivatives of the spherical Bessel functions are now implemented as well (and accessible using the derivative optional argument). This completes the implementation of #3113.

[rgommers]

Looks promising. You may want to mention this on the mailing list - there may be people interested in trying it out.

I see in gh-3113 that @certik originally asked for this, so pinging him here.

[certik]

Thanks for working on it @tpudlik. I looked at the tests and I think that's exactly what I wanted!

[argriffing]

This looks useful, but I'm not going to review the technical details. The PR has some python comments related to technical limitations regarding 'complex infinity', and I think it is handled well. I hope this can be merged soon instead of waiting indefinitely for a reviewer to dive into DLMF.

The performance improvement is not dramatic, at a little over an order of magnitude

That sounds more than dramatic enough!

[ev-br]

(I did not review the details). I wonder if we should

• test these against boost (using the same tests as for the sph_* functions)
• make sure these pass the tests for sph_*.
• deprecate the sph_* functions in favor of these.

[ev-br]

ev-br@e305b53 tries to add a quick-and-dirty testing of spherical_jn and spherical_kn against boost. Both fail :-(. OTOH, sph_jn fails as well.

[tpudlik]

I took a look at the boost tests---thank you for bringing them to my attention, @ev-br! I was able to get them to pass by fixing a couple bugs that reduced accuracy and raising the tolerance to 1e-13 for spherical_jn.

I also added the tests for sph_* from test_basic.py.

[codecov-io]

@@            master   #5623   diff @@
======================================
  Files          237     238     +1
  Stmts        43523   43542    +19
  Branches      8162    8166     +4
  Methods          0       0       
======================================
+ Hit          33943   33962    +19
  Partial       2601    2601       
  Missed        6979    6979       

Review entire Coverage Diff as of 372b9cc

Powered by Codecov. Updated on successful CI builds.

[ev-br]

Three more general comments (hope you don't mind piecewise review too much):

• Testing against boost seems useful for Jn and Yn. Is there a way of testing spherical_in and spherical_kn in a similar way: Boost, GSL, Mathematica? It's possible that there is no readily available alternative implementation, so maybe spend a limit amount of effort on this.
• What do we do with the existing sph_* functions? One option is to deprecate them in favor of the new ones; another one is to just remove the old ones and rename spherical_* to sph_*.
• I'm two minds on having a boolean derivative keyword. Maybe better make it integer, so that higher derivatives can be implemented later.

I'm +0 on all these points above, so that none of these is blocking.

[pv]

These functions can probably be tested against mpmath, cf. tests in test_mpmath.py.

[tpudlik]

I deprecated the sph_* functions. Simply replacing them with the new ones seems undesirable because the change would not be backwards-compatible. The new and old functions' signatures are different: sph_* returns a (2, n + 1) array for scalar inputs, which includes both the function values and the derivatives for all orders up to n, while the new functions exhibit standard ufunc behavior (scalar inputs give a scalar output).

The derivative boolean keyword can be replaced with an integer counter in a backwards-compatible way later if higher derivatives become available, since Python casts False to a 0 and True to a 1. (In fact, whether we use integers of booleans now is a feature of the documentation, not the implementation, which works just as well if passed 0 or 1.) But the introduction of higher derivatives is unlikely: to the best of my knowledge, the recurrence relations that allow the direct computation of the first derivatives (http://dlmf.nist.gov/10.51.E2 and E5) don't have analogs for higher ones.

As far as testing is concerned: boost and Mathematica unfortunately do not include the modified spherical Bessel functions. But I added tests to test_mpmath.py. Like for the cylindrical Bessel functions, the domain of the argument must sometimes be restricted (the functions are oscillatory, and high precision cannot be achieved for very large argument) and I restrict the order of the functions to 200 (150 for spherical_kn for real arguments, which sometimes calls overflow a little too early for larger order). Within these constraints all functions pass the tests, except for spherical_kn for a few complex arguments with much smaller imaginary than real parts, where a loss of precision occurs. I marked the latter as a "known failure".

[pv]

Regarding the failure of complex spherical_kn, I think the correct procedure here is

1. try to fix it
2. if this is too hard, and you know where on the complex plane the failure is, the function should return nan and emit (sf_error.NO_RESULT, "total loss of precision"), instead of returning known-incorrect results.

[tpudlik]

Just a quick update: I'm working on spherical_kn. I have some ideas for fixing it, and if that doesn't work, I'll add warnings in the relevant domain. I'm hoping to have something to show by next week.

[tpudlik]

I've looked into the spherical_kn problem more closely. I have been able to find only one region in which there is loss of precision: for n = 1, when z = -1 + 1j*e for e small (less than 1e-5). The cause of the problem is a silent precision loss in kv. For example,

>>> scipy.special.kv(1.5, -1 - 1j*1e-20)
(-5.6498675426425281e-17+2.2204460492503131e-16j)
>>> mpmath.besselk(1.5, -1 - 1j*1e-20)
mpc(real='-3.4068610448155485e-20', imag='-1.7034305224077743e-40')

I could add an error to spherical_kn, but I think it makes more sense to do so in kv. I worry that there are more such difficult points in the complex plane (and for other orders) and I just haven't found them yet. For now I created an issue, #5844.

[pv]

Probably not needed? The spherical functions are included in the _ufunc.so extension?

[pv]

Ok, in that case probably not so much can be done in this PR to address the K issue.
Final technical comments: spherical_bessel.py -> _spherical_bessel.py to mark it private.
Otherwise LGTM.

[tpudlik]

We should be all set now! (Except Travis CI doesn't seem to have started a build on my commit---is there a way to trigger it manually?)

[rgommers]

@tpudlik the PR needs to be rebased. I'm not sure if TravisCI tests will run then (I've never seen them not trigger at all), but otherwise maybe the simplest solution is to close this PR and open a new one with the same content. I'm not aware of a way to manually trigger.

[larsoner]

I think Travis won't trigger if the PR can't be merged. So a rebase should fix it.

[pv]

thanks, merged