Compatibility with QuadGK

[aurelio-amerio]

Hello, the update of QuadGK (a popular library for numerical integration) from v2.5 to v2.6 broke compatibility with Unitful.

To restore compatibility, one would have to extend the QuadGK function cachedrule to also work with quantities with something like

@generated function cachedrule(::Type{T}, n::Integer) where {T<:AbstractQuantity}
    TF = typeof(float(real(one(T))))
    cache = haskey(rulecache, TF) ? rulecache[TF] : (rulecache[TF] = Dict{Int,NTuple{3,Vector{TF}}}())
    :(haskey($cache, n) ? $cache[n] : ($cache[n] = kronrod($T, n)))
end

Unitful and QuadGK don't have each other as dependencies, so it would not be possible to extend either without adding more dependencies.

Do you think it would be beneficial to create a new julia package (like UnitfulQuad or something similar) to make Unitful talk nicely with QuadGK, or do you have any suggestions on how to proceed? I think out of the box numerical integration with units is a strong selling point for Julia!

[sostock]

There is already a package that is supposed to handle compatibility between Unitful and QuadGK. It is called UnitfulIntegration. However, it hasn’t been updated in a while and is not compatible with Unitful 1.x. My suggestion would be to bring UnitfulIntegration up to date to support the current versions of QuadGK and Unitful.

For Julia ≥ 1.10 1.9, we could add the necessary methods as a package extension (JuliaLang/julia#47695). UnitfulIntegration would then no longer be needed.

[aurelio-amerio]

Alright! I'm working on updating the package and I will issue a pull request as soon as possible!

[aurelio-amerio]

I have submitted a pull request on UnitfulIntegration

[stevengj]

Can you give a specific test case?

I just tried it with the latest versions (QuadGK 2.6.0 and Unitful v1.12.2) in Julia 1.8.2, and it works fine for me:

julia> using Unitful, QuadGK

julia> f(t) = 1.0u"J"
f (generic function with 1 method)

julia> quadgk(f, 0.0u"s", 10.0u"s")
(10.0 J s, 0.0 J s)

[stevengj]

It seems like there is a bug in Julia where if you do using Unitful, QuadGK it works fine, but if you do using QuadGK, Unitful it fails: JuliaLang/julia#48295

[vtjnash]

there are a couple things here that are forbidden to generated functions.

@generated function cachedrule(::Type{T}, n::Integer) where {T<:AbstractQuantity}
    TF = typeof(float(real(one(T))))
    cache = haskey(rulecache, TF) ? rulecache[TF] : (rulecache[TF] = Dict{Int,NTuple{3,Vector{TF}}}())
    :(haskey($cache, n) ? $cache[n] : ($cache[n] = kronrod($T, n)))
end

1. calling one is UB, since it is expected that the user may overload it with new methods
2. accessing mutable state (e.g. rulecache here) is UB, as it may run simultaneously on multiple threads or on no threads at all in this process

For addressing 1, the solution is simply to move that code outside of the generated function, so that the argument to the generated function is TF instead of T

[sostock]

Infinite bounds do not work, even with using Unitful, QuadGK:

julia> using Unitful, QuadGK

julia> quadgk(x->exp(-x/(1.0u"m")), 0.0u"m", Inf*u"m", atol=0.0u"m")
ERROR: DomainError with Inf m:
integrand produced NaN m in the interval (0.0 m, Inf m)
Stacktrace:
[...]

[stevengj]

Update: this turned out to be a bug in the use of @generated in QuadGK, should be closed by JuliaMath/QuadGK.jl@1959d34 (I'll tag a release soon), so this issue can be closed I think?

[stevengj]

Closed by JuliaRegistries/General#75775

[sostock]

Infinite bounds still don’t work:

(jl_s3h7TF) pkg> st
Status `/tmp/jl_s3h7TF/Project.toml`
  [1fd47b50] QuadGK v2.7.0 `https://github.com/JuliaMath/QuadGK.jl#master`
  [1986cc42] Unitful v1.12.2

julia> using Unitful, QuadGK

julia> quadgk(x->exp(-x/(1.0u"m")), 0.0u"m", Inf*u"m")
ERROR: DomainError with Inf m:
integrand produced NaN m in the interval (0.0 m, Inf m)
Stacktrace:
[...]

[stevengj]

The problem with infinite bounds is a completely unrelated problem. I would suggest closing this issue in favor of a more specialized issue JuliaMath/QuadGK.jl#64

[stevengj]

Should be fixed now, in any case. On master (soon to be version 2.7.0):

julia> quadgk(x->exp(-x/(1.0u"m")), 0.0u"m", Inf*u"m")
(1.0 m, 4.507383379289404e-11 m)

[sostock]

Great, thanks!