JuliaMath / QuadGK.jl

adaptive 1d numerical Gauss–Kronrod integration in Julia
MIT License
272 stars 37 forks source link

infinite limits with units and segbuf broken #95

Open lxvm opened 11 months ago

lxvm commented 11 months ago

The issue is that segbufs don't always work with infinite limits, such as when the limits have units. This is a MWE

using QuadGK, Unitful

Tu = typeof(1.0u"m")
segbuf = QuadGK.alloc_segbuf(Tu, Tu, Tu)

quadgk(x -> 1/(1+oneunit(x)/x), 0.0u"m", 1.0u"m", segbuf=segbuf) # OK (0.3068528194400547 m, 1.9931833961095435e-11 m)

quadgk(x -> 1/(1+oneunit(x)/x), 0.0u"m", Inf*u"m", segbuf=segbuf) # broken

which gives the error

ERROR: DimensionError: m and 0.0 are not dimensionally compatible.
Stacktrace:
  [1] convert(#unused#::Type{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, x::Float64)
    @ Unitful ~/.julia/packages/Unitful/jiUeO/src/conversion.jl:106
  [2] QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}(a::Float64, b::Float64, I::Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, E::Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}})
    @ QuadGK ~/.julia/dev/QuadGK/src/evalrule.jl:3
  [3] convert(#unused#::Type{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, s::QuadGK.Segment{Float64, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}})
    @ QuadGK ~/.julia/dev/QuadGK/src/evalrule.jl:10
  [4] setindex!(A::Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}, x::QuadGK.Segment{Float64, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, i1::Int64)
    @ Base ./array.jl:969
  [5] macro expansion
    @ ./broadcast.jl:974 [inlined]
  [6] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [7] copyto!
    @ ./broadcast.jl:973 [inlined]
  [8] copyto!
    @ ./broadcast.jl:926 [inlined]
  [9] materialize!
    @ ./broadcast.jl:884 [inlined]
 [10] materialize!
    @ ./broadcast.jl:881 [inlined]
 [11] do_quadgk(f::QuadGK.var"#18#27"{var"#11#12", Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}, s::Tuple{Float64, Float64}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}})
    @ QuadGK ~/.julia/dev/QuadGK/src/adapt.jl:39
 [12] #50
    @ ~/.julia/dev/QuadGK/src/adapt.jl:235 [inlined]
 [13] handle_infinities(workfunc::QuadGK.var"#50#51"{Nothing, Nothing, Int64, Int64, typeof(LinearAlgebra.norm), Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}}}, f::var"#11#12", s::Tuple{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}})
    @ QuadGK ~/.julia/dev/QuadGK/src/adapt.jl:126
 [14] quadgk(::Function, ::Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, ::Vararg{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}; atol::Nothing, rtol::Nothing, maxevals::Int64, order::Int64, norm::Function, segbuf::Vector{QuadGK.Segment{Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, Quantity{Float64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}}})
    @ QuadGK ~/.julia/dev/QuadGK/src/adapt.jl:234
 [15] top-level scope
    @ REPL[11]:1