Comparing besselj speed and values between SpecialFunctions.jl and Bessels.jl

BesseljComparison.jl

using SpecialFunctions, Bessels, LinearAlgebra, Plots

const spj = SpecialFunctions.besselj
const bej = Bessels.besselj

function foo!(A, B, f::T, ns, ms, N) where T
  for (j, n) in enumerate(ns)
    for (i, m) in enumerate(ms)
      A[i, j] = 0
      t = @elapsed for _ in 1:N
        A[i, j] = max(A[i, j], abs(f(n, 10.0^(m + rand()))))
      end
      B[i, j] = t / N
    end
  end
end
const ns = -50:50
const ms = -64:31
const A1 = zeros(length(ms), length(ns))
const A2 = zeros(length(ms), length(ns))
const A3 = zeros(length(ms), length(ns))
const B1 = zeros(length(ms), length(ns))
const B2 = zeros(length(ms), length(ns))
const B3 = zeros(length(ms), length(ns))

foo!(A1, B1, spj, ns, ms, 1)
foo!(A2, B2, bej, ns, ms, 1)
foo!(A3, B3, (n, x) -> (spj(n, x) - bej(n, x)), ns, ms, 1)
const M = 1000
t1 = @elapsed foo!(A1, B1, spj, ns, ms, M)
t2 = @elapsed foo!(A2, B2, bej, ns, ms, M)

function bar(n, x)
  a = spj(n, x)
  b = bej(n, x)
  return (a - b) / max(abs(a), abs(b))
end
t3 = @elapsed foo!(A3, B3, bar, ns, ms, M)

h1 = heatmap(ns, ms, log10.(abs.(A1)), title="SpecialFunctions")
ylabel!("Exponent")
h2 = heatmap(ns, ms, log10.(abs.(A2)), title="Bessels")
h3 = heatmap(ns, ms, log10.(abs.(A3)), title="Normalised Difference")
xlabel!("Order")
ylabel!("Exponent")

indices = findall(-3 .<= ns .<= 3)
h4 = plot(bg=:white, legend=:topleft)
for i in indices
  plot!(h4, ms, log10.(abs.(A3[:, i]))[:], label="$(ns[i])")
end
ylabel!("Log10 |Δ|")
xlabel!("Exponent")
plot(h1, h2, h3, h4, layout=@layout [a b; c d])
savefig("BesseljComparison.pdf")

h1 = heatmap(ns, ms, log10.(B1), title="SpecialFunctions")
xlabel!("Order")
ylabel!("Exponent")
h2 = heatmap(ns, ms, log10.(B2), title="Bessels")
xlabel!("Order")
ylabel!("Exponent")
h3 = heatmap(ns, ms, log10.(B2 ./ B1), title="Bessels / SpecialFunctions")
xlabel!("Order")
ylabel!("Exponent")

y = (B2[1, :] ./ B1[1, :])[:]
(_, ind) = findmax(y)
x = ns[ind]
h4 = plot(bg=:white, legend=:topleft)
plot!(h4, ns, log10.(y), label=nothing)
plot!(h4, [x, x], log10.(collect(extrema(y))), label="$x")
ylabel!("Time ratio")
xlabel!("Order")
title!("Exponent = $(ms[1])")

plot(h1, h2, h3, h4, layout=@layout [a b; c d])
savefig("BesseljTime.pdf")