switch to Bessels.jl for real bessel functions #409
Comments
|
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 Of course, I think Bessel functions are special enough to include here, but there would be a couple of concerns I would have.
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. |
|
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. |
|
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.
Are there any benchmarks one can check?
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. |
This seems like a really great model to follow. I would be very interested in something like that to happen here. Particular for
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 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.
I agree with this! I have opened an issue at JuliaMath/Bessels.jl#41. This really only affects two functions ( |
|
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: This is some amazing work @heltonmc ! |
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 JuliaDiff/DiffRules.jl#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 #160 (comment)) 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. |
I agree. That is the best approach when
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 |
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. |
|
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 #409 (comment) and #409 (comment) is better. |
|
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. |
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. |
|
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). |
|
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. |
Here is a first shot for an automated way ( # 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 This could be repeated for |
|
Yes, |
|
We could use git subtrees or git submodules for a built-in system to handle this. |
|
Why do you need subtrees/submodules at all if the code is being moved to another package? |
|
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. |
|
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. |
|
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. |
|
I was talking about compilation latency after loading the package, not precompilation before even using it. |
|
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. |
|
Latency seems quite good. I wonder if we should just replace the ccalls to amos with calls to Bessels.jl, the topic of this issue. I was asking @oscardssmith about doing this and he pointed me to this issue. Reading the discussion here, I feel that doing so would generally be a good idea. A lot of this goes back to primordial Julia, but today, we would actually have an |
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.
The text was updated successfully, but these errors were encountered: