switch to Bessels.jl for real bessel functions

[oscardssmith]

https://github.com/JuliaMath/Bessels.jl has pure Julia bessel functions for real arguments and orders that are frequently significantly faster than the ones currently used (provided by AMOS). The package was just registered, but the implementations are solid and well tested. Also, on the roadmap for the future are the complex orders and arguments.

[heltonmc]

Thanks for starting this discussion! I am very interested to see everyone's thoughts on this and open to whatever the consensus of the community is on this and would volunteer time to make this happen.

Though, it does appear that discussions like this have frequently occurred here and on discourse about what functions to include in SpecialFunctions.jl, the future of the package, finding maintainers, and reviewing PRs. It does appear like there is a lot of uncertainty around PRs about if they should be included or just their general quality.

Of course, I think Bessel functions are special enough to include here, but there would be a couple of concerns I would have.

1. Testing. Currently Bessels.jl tests many values against the AMOS routines provided by SpecialFunctions.jl. We would have to figure out a different way to test (ideally still having the AMOS routines available in a different package - openspeclib?).
2. Many packages depend on SpecialFunctions.jl which could create some instability for a large part of the julia ecosystem. The code in Bessels.jl to my knowledge is much faster (3x-10x) but it does take a lot longer to compile, is still new and being tested, and being refactored. Inclusion of Bessels.jl would probably over double the source code and increase maintainer burden.
3. Intermixing native julia code and non julia code. One of the advantages I think of Bessels.jl is it has no dependencies and is written entirely in Julia.

All of these can of course be alleviated with some effort, but it may also be worth considering a different approach. I think (speaking very idealistic) it would be great if SpecialFunctions.jl was able to provide only native julia implementations whereas we could put the fortran libraries into some other package (openspeclib.jl?) to test against.

However, it might be best to have SpecialFunctions.jl serve as an umbrella package that could reexport smaller lighter-weight dependencies. For example, I'm sure lots of people might just need a lightweight gamma functions but don't need Bessel or elliptic functions. This also helps with CI and streamlining testing. I think SpecialFunctions.jl could serve as a home and major resource to point users to these packages that were under some sort of purveyorship of the JuliaMath group. A great example is (https://juliadiff.org/) that provides a great overview of the automatic differentiation tools in the ecosystem.

I would then pose that classes of special functions be their own separate packages that consists of groups of similar functions: 1. gamma, inverse gamme, incomplete gamma, polygamma, zeta, 2. bessel, airy, hankel, spherical bessel functions, 3. kelvin, lommel, Anger, Weber, Struve, Scorer, 4. hypergeometric functions, 5. elliptic function, 6. orthogonal polynomials.

Now several of these packages already exist (e.g. HypergeomtricFunctions.jl, Elliptic.jl, ClassicalOrthogonalPolynomials.jl, Bessels.jl) but it would be nice if SpecialFunctions.jl served as a more centralized point and curator of the julia mathematical ecosystem. Of course this takes a lot of work, but I would be happy to volunteer time right now to make that happen.

Edited for conciseness.

[stevengj]

I think it would make a lot of sense to have some subsidiary packages for particular special functions, while keeping this package as an umbrella package (but changing it to simply import from the specialized packages), and encourage downstream dependencies to require the subsidiary package in order to lighten their dependencies.

[devmotion]

I agree, I think moving some classes of special functions to (or keeping them in) separate packages and only reexporting them in SpecialFunctions makes sense. We moved the log/exp functions from StatsFuns to a separate package LogExpFunctions a while ago and it seems development became more active, but the biggest benefit was probably that it made it possible for many packages to depend on and use LogExpFunctions that would never have taken a dependency on StatsFuns (and Rmath).

Generally, I guess it would be simpler though to extract things to separate packages before refactoring or rewriting them (e.g., in native Julia code). That would allow to preserve git history (we did that in LogExpFunctions) and would ensure that no breaking changes are accidentally introduced in SpecialFunctions during the split.

The code in Bessels.jl to my knowledge is much faster (3x-10x)

Are there any benchmarks one can check?

but it does take a lot longer to compile

Can we fix the compile time issues before using Bessels.jl in SpecialFunctions? SpecialFunctions is a central package in the ecosystem (currently 2181 direct and indirect dependents), so it would be very unfortunate to introduce any major compile time or runtime regressions.

[heltonmc]

We moved the log/exp functions from StatsFuns to a separate package LogExpFunctions a while ago and it seems development became more active

This seems like a really great model to follow. I would be very interested in something like that to happen here. Particular for gamma and related functions as it looks like a lot of those already have julia implementations. We could then port the Bessels.gamma(x::Float64) implementation over to there and have Bessels.jl depend on that instead.

Are there any benchmarks one can check?

Of course! I used that qualifier because performance really depends on the specific argument and which branch you fall in. Take for example the Hankel function (which has a lot of different branches in both AMOS and Bessels.jl)

julia> @benchmark Bessels.hankelh1(2.3, x) setup=(x=rand()*100)
BenchmarkTools.Trial: 10000 samples with 940 evaluations.
 Range (min … max):  102.350 ns … 805.496 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     103.856 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   178.015 ns ± 182.483 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █    ▂▃▂               ▁                         ▁▂  ▁ ▁ ▁▁▁  ▁
  █▇▅▅████▄▁▁▃▁▁▁▁▁▁▁▁▁▁▁█▅▄▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁███████████████ █
  102 ns        Histogram: log(frequency) by time        731 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark SpecialFunctions.hankelh1(2.3, x) setup=(x=rand()*100)
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max):  1.175 μs …   4.062 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.404 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.581 μs ± 366.964 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

        ▁██▆▅                                                  
  ▂▂▂▂▃▆██████▇▅▄▄▄▄▄▃▂▂▂▁▁▁▂▂▂▁▁▁▂▁▁▂▁▂▂▂▂▂▂▂▂▂▃▆▇▆▅▄▃▂▂▂▂▂▂ ▃
  1.18 μs         Histogram: frequency by time         2.5 μs <

 Memory estimate: 32 bytes, allocs estimate: 2.

Of course the longest branch is still shorter than the shortest branch in the amos routine, but it just depends. Looking at an example that doesn't allocate (MPFR routine) for very large orders..

julia> @benchmark SpecialFunctions.bessely(1e6, x) setup=(x=rand()*1e4 + 1e6 )
BenchmarkTools.Trial: 2254 samples with 1 evaluation.
 Range (min … max):  2.189 ms …  2.452 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.222 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.216 ms ± 20.992 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▂▇    ▂                   ▅ █                              
  ▆██▃▃▅▅█▅▃▄▃▃▃▃▃▃▃▃▃▃▃▅▅▆▅▄█▇██▄▃▄▄▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂ ▃
  2.19 ms        Histogram: frequency by time        2.27 ms <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark Bessels.bessely(1e6, x) setup=(x=rand()*1e4 + 1e6 )
BenchmarkTools.Trial: 10000 samples with 909 evaluations.
 Range (min … max):  119.637 ns …   2.546 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     120.783 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   238.958 ns ± 394.606 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █▆▃                                                           ▁
  ███▄▁▁▁▁▆▇▆▆▆▆▆▆▇▅▆▅▆▆▆▇▅▆▆▇▆▆▆▅▆▆▆▆▆▆▇▆▆▆▆▆▆▆▇▆▆▆▆▅▆▅▆▇▆▆▆▆▆ █
  120 ns        Histogram: log(frequency) by time       2.28 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

So again there is potential for a large amount of performance gains depending on the use case. Maybe a more widely used example with only 2 branches would be besselk of order 0.

julia> @benchmark Bessels.besselk(0, x) setup=(x=rand()*100)
BenchmarkTools.Trial: 10000 samples with 997 evaluations.
 Range (min … max):  19.642 ns … 43.214 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     26.119 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   26.119 ns ±  0.932 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

      ▁                                  ▁▆█▄                 ▁
  ▅▄▃▁█▆▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄████▅▄▅▆▅▅▄▃▄▃▅▃▁▃▄▇ █
  19.6 ns      Histogram: log(frequency) by time      28.9 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

julia> @benchmark SpecialFunctions.besselk(0, x) setup=(x=rand()*100)
BenchmarkTools.Trial: 10000 samples with 723 evaluations.
 Range (min … max):  173.351 ns …  2.379 μs  ┊ GC (min … max): 0.00% … 89.15%
 Time  (median):     228.043 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   250.504 ns ± 63.841 ns  ┊ GC (mean ± σ):  0.17% ±  1.27%

        ▄█▄                                                     
  ▁▁▁▁▁▁████▇▆▄▃▂▂▂▄▄▄▃▃▂▂▂▂▂▂▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  173 ns          Histogram: frequency by time          502 ns <

 Memory estimate: 16 bytes, allocs estimate: 1.

So there are definitely a lot of potential advantages in runtime for these as well as none of the function allocate.

Can we fix the compile time issues before using Bessels.jl in SpecialFunctions?

I agree with this! I have opened an issue at https://github.com/JuliaMath/Bessels.jl/issues/41. This really only affects two functions (besselj and bessely) and is limited to one function that completely unrolls an asymptotic expansion. The first run for these functions is definitely noticeable (2s) but I will work on this today to improve. Edit: the culprit functions accounts for over half the first run time so hopefully we can find a way to improve this so first run is under 1 second.

[t-bltg]

Cc @juliohm for https://github.com/JuliaPlots/UnicodePlots.jl/issues/291.

[juliohm]

An alternative path would be to develop a separate package named SpecialFuns.jl (or any similar name), which supersedes Bessels.jl and related functionality without dependencies or any promise of backwards compatibility. That way maintenance would be easier for the people working on it, and the switch from SepcialFunctions.jl to SpecialFuns.jl would be trivial downstream. No playing with version numbers, but a simple string replacement s/SpecialFunctions/SpecialFuns/g

Unless SpecialFunctions.jl aims to be a lightweight (zero-dependency) package for special functions, the effort of plugging Bessels..jl into it doesn't seem to pay off in the long-term. Of course I may be wrong.

[ParadaCarleton]

An alternative path would be to develop a separate package named SpecialFuns.jl (or any similar name), which supersedes Bessels.jl and related functionality without dependencies or any promise of backwards compatibility. That way maintenance would be easier for the people working on it, and the switch from SepcialFunctions.jl to SpecialFuns.jl would be trivial downstream. No playing with version numbers, but a simple string replacement s/SpecialFunctions/SpecialFuns/g

Unless SpecialFunctions.jl aims to be a lightweight (zero-dependency) package for special functions, the effort of plugging Bessels..jl into it doesn't seem to pay off in the long-term. Of course I may be wrong.

There would be a huge payoff. Right now, automatic differentiation doesn't work with Bessel functions, which is a huge problem. There's about half a dozen distributions that HMC can't work with in Turing alone. Switching to Bessels.jl would be a huge improvement.

That being said, it looks like compilation time has already been fixed, amazingly enough:

julia> @time using Bessels
  0.021716 seconds (37.68 k allocations: 3.757 MiB)

julia> @time besselix(1, 1)
  0.000001 seconds

This is some amazing work @heltonmc !

[devmotion]

automatic differentiation doesn't work with Bessel functions

SpecialFunctions has contained ChainRules definitions for Bessel functions for quite some time, and some bug fixes in ForwardDiff and ReverseDiff made it possible to finally merge https://github.com/JuliaDiff/DiffRules.jl/pull/79 yesterday which added similar functionality to DiffRules-based AD backends such as ForwardDiff and ReverseDiff. Of course, that won't help if you want to differentiate Bessel functions wrt to the order (see eg https://github.com/JuliaMath/SpecialFunctions.jl/issues/160#issuecomment-755375109) but basically everything else should work now.

Generally, I also assume the best approach for AD in this case are custom rules, regardless of whether or not the primal function is implemented in native Julia or not.

[heltonmc]

Generally, I also assume the best approach for AD in this case are custom rules, regardless of whether or not the primal function is implemented in native Julia or not.

I agree. That is the best approach when nu=0 as they can be expressed with a single call for nu=1 which is also optimized. I originally thought the oscillating functions besselj and bessely might be less accurate for arbitrary orders due to the subtraction. But after some quick tests it looks like they can be calculated around ~5 ULP error when the Bessel functions are correctly rounded. There was a few cases where I saw Ulp > 40..... One thing to consider is the derivative then requires 2 Bessel function calls. If ForwardDiff could potentially be just as accurate and give the derivative in one call that could be an advantage.

Of course, that won't help if you want to differentiate Bessel functions wrt to the order (see eg #160 (comment)) but basically everything else should work now.

Also, just to mention that Bessels.jl doesn't even allow for dual numbers right now so doesn't work with ForwardDiff. This should be a simple fix to just allow both Float64 and Dual types to propagate. I don't have that much experience with AD so that would be a very nice addition to Bessels.jl if anyone would like to contribute. I should say that we shouldn't also just expect to drop in AD and get accurate derivatives wrt to order. I think besselj and besseli would be the easiest as their power series works for integer and non-integer orders. And I just checked that in the case of besselj the derivative also has stable recurrence... these two functions would be interesting cases to test. I'm less confident for the modified forms (https://github.com/cgeoga/BesselK.jl) as we depend on different algorithms for integer nu and non-integer nu. I'm unsure if derivative information around those branches would be conserved. Though, this would be a very interesting investigation and a great addition!

[heltonmc]

An alternative path would be to develop a separate package named SpecialFuns.jl (or any similar name), which supersedes Bessels.jl and related functionality without dependencies or any promise of backwards compatibility.

I do particularly like this idea.... I think the challenge would be finding maintainers. It would make sense if @devmotion and some of the contributors here could help maintain. But they are also some of the most prolific contributors in the entire Julia ecosystem. I'm not really sure how much time they would have to take on additional things considering it would be a big project. I would like to see the functions here in gamma.jl be pruned into a new repository with no external dependencies. I briefly played around with that trying to maintain git history but couldn't quite get it right. If someone wants to do that I would be happy to maintain it. Bessels.jl has it's own gamma implementation for Float64 that is accurate to <1.5 ULP and I'm working on loggamma for real arguments with that same accuracy. That package would have a lot of nice functionality already that the proposed new special functions could be built on.

[devmotion]

I don't think a separate SpecialFunctions package is the best way to go - I assume since SpecialFunctions is such a core part of the Julia ecosystem in most cases downstream projects and packages will only end up depending on both SpecialFunctions and SpecialFunctions2, im- or explicitly. And if there's a clearly better implementation, why should we keep the other one around instead of letting the whole Julia ecosystem benefit from the improved version? I think the approach outlined in https://github.com/JuliaMath/SpecialFunctions.jl/issues/409#issuecomment-1212641119 and https://github.com/JuliaMath/SpecialFunctions.jl/issues/409#issuecomment-1212844459 is better.

[devmotion]

And yeah, it's already hard finding maintainers and reviewers for SpecialFunctions, so I assume it would become even more difficult if some focus on one, and some on the other alternative of SpecialFunctions.

[heltonmc]

I think the approach outline in #409 (comment) and #409 (comment) is better.

Yes, sorry I wasn't clear I would be very supportive of that outline and would volunteer to help maintain here and other packages. One thing is figuring out tests as we test against SpecialFunctions (I enjoy having the Amos routines so easily accessible but I'm sure I can set up my tests directly to openspecfun). I would be in favor though of separating out Julia code and having SpecialFunctions.jl bridge the native Julia and external libraries until eventually those are replaced. Ideally, Bessels.jl could depend on a lightweight hypothetical Gamma.jl package that SpecialFunctions could also reexport.

[stevengj]

Onet path forward would be to split more functionality out of SpecialFunctions.jl into subpackages, e.g. GammaFunctions.jl, ErrorFunctions.jl, ExponentialIntegrals.jl, which are loaded and re-exported by SpecialFunctions.jl for backward compatibility. (i.e. the total amount of code wouldn't change).

Then we could encourage downstream packages to rely on only the subpackage(s) that they need rather than on all of SpecialFunctions, and add new special functions in the future to their own subpackage without necessarily including them in SpecialFunctions.jl (so that this package stops growing larger and larger).

[devmotion]

To reemphasize (I think I mentioned it above), I think it would be very helpful if the specialized subpackages would preserve the git history of the files extracted from SpecialFunctions. I remember quite a few relatively recent PRs that fixed subtle bugs that only became apparent when we switched more parts of StatsFuns to native Julia implementations based on SpecialFunctions. It would be unfortunate if the reason and motivation for these fixes and deviations from papers/fortran code/etc would be lost in the transition.

[t-bltg]

I briefly played around with that trying to maintain git history but couldn't quite get it right.

I think it would be very helpful if the specialized subpackages would preserve the git history of the files extracted from SpecialFunctions

Here is a first shot for an automated way (bash on linux) to extract a set of files from SpecialFunctions while preserving the git history, to a new repo named https://github.com/JuliaMath/GammaFunctions.jl.

# sudo apt -y install git-filter-repo

pkg_name=GammaFunctions  # new package name

# extracted files from `SpecialFunctions.jl`
keep+=(src/gamma.jl)
keep+=(test/gamma.jl)
keep+=(test/gamma_inc.jl)

check() {
  # check that some random commits affecting retained files are found (history)
  # TODO: more checks ?
  git log | grep 'Use `mpfr_lngamma` for `loggamma(::BigFloat)`'
  git log | grep 'make hurwitz zeta(s,z) more Julian: require complex inputs to get complex outputs'
  git log | grep 'split tests into components'
}

#################################################
special_functions=$(mktemp -d)  # temporary clone

git clone https://github.com/JuliaMath/SpecialFunctions.jl $special_functions

# extract commits affecting kept files
(
  cd $special_functions

  git filter-repo ${keep[@]/#/--path }
  echo "extracted $(git rev-list --count master) commit(s)"
  check
)

# create julia package structure
julia -e '
  using PkgTemplates
  tpl = Template(
    user = "JuliaMath",
    dir = pwd(),
    julia = v"1.3.0",
    plugins = [
      GitHubActions(),
      Codecov(),
      # TODO: add more plugins ?
    ]
  )
  tpl(first(ARGS))
' $pkg_name

# merge selected commits into new package
(
  cd $pkg_name

  git remote add temp $special_functions
  git fetch temp master
  git merge -m 'merge selected files' --allow-unrelated-histories FETCH_HEAD
  git remote remove temp
  check
  git remote -v
)

I've used this procedure with git-filter-repo to move RecipesBase and RecipesPipeline into Plots for consolidation into a monorepo in https://github.com/JuliaPlots/Plots.jl/pull/4419, while rewriting Plots history to include the history of these modules, and we've had no issue with that approach.

This could be repeated for ErrorFunctions.jl, ExponentialIntegrals.jl, ...

[stevengj]

Yes, git filter-repo is the way to go here. This can prune out the files at a course level, and any remaining duplication can be edited out manually.

[ParadaCarleton]

We could use git subtrees or git submodules for a built-in system to handle this.

[giordano]

Why do you need subtrees/submodules at all if the code is being moved to another package?

[ParadaCarleton]

I'm actually not sure if separating this into separate packages is even necessary post-1.9, given the big improvements in compile times. Especially given there's a lot more on the way already, and almost all workflows in Julia are going to require loading quite a bit of the code here anyways.

[giordano]

My understanding was that the separation was suggested in order to ease maintenance of code and speed up development. I don't see how splitting up the code into multiple packages would help with compilation latency even before v1.9.

[oscardssmith]

The reason splitting speeds up latency is that it means that packages that only rely on a subset of the functionality don't have to load all the code. This benefit is enhanced with 1.9 since in 1.9 multiple packages will increase the efficiency of parallel precompilation.

[giordano]

I was talking about compilation latency after loading the package, not precompilation before even using it.

[oscardssmith]

The same thing applies for packages that only require some of the functions in terms of load speed. Loading more code takes more time and memory. It's a lot faster to not load as much code if you don't need to.