Last active
October 17, 2022 20:00
-
-
Save jwscook/4712f7120870f32378f00e69781e604a to your computer and use it in GitHub Desktop.
Comparing besselj speed and values between SpecialFunctions.jl and Bessels.jl
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment