JuliaPhysics/Measurements.jl

Nested Types issue

DanBoland opened this issue · 3 comments

DanBoland commented

I'm having trouble integrating with units and measurements.

using Unitful, QuadGK, Measurements
f(x)=x^2
quadgk(f, 0u"m", (10±1)u"m") # fails: QuadGK + Measurements + Unitful

that fails with:

MethodError: no method matching Float64(::Measurements.Measurement{Float64})

convert(::Type{Float64}, ::Measurements.Measurement{Float64})@number.jl:7
setindex!(::Array{Float64,1}, ::Measurements.Measurement{Float64}, ::Int64)@array.jl:847
eignewt(::Array{Measurements.Measurement{Float64},1}, ::Int64, ::Int64)@gausskronrod.jl:43
kronrod(::Type{Measurements.Measurement{Float64}}, ::Int64)@gausskronrod.jl:193
macro expansion@gausskronrod.jl:257[inlined]
cachedrule@gausskronrod.jl:257[inlined]
cachedrule@gausskronrod.jl:249[inlined]
do_quadgk(::Main.workspace192.var"#1#2", ::Tuple{Unitful.Quantity{Measurements.Measurement{Float64},𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},Unitful.Quantity{Measurements.Measurement{Float64},𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}}}, ::Int64, ::Nothing, ::Nothing, ::Int64, ::typeof(LinearAlgebra.norm))@adapt.jl:7
(::QuadGK.var"#28#29"{Nothing,Nothing,Int64,Int64,typeof(LinearAlgebra.norm)})(::Function, ::Tuple{Unitful.Quantity{Measurements.Measurement{Float64},𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},Unitful.Quantity{Measurements.Measurement{Float64},𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}}}, ::Function)@adapt.jl:179
handle_infinities@adapt.jl:113[inlined]
#quadgk#27@adapt.jl:177[inlined]
quadgk(::Function, ::Unitful.Quantity{Measurements.Measurement{Float64},𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}}, ::Unitful.Quantity{Measurements.Measurement{Float64},𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}})@adapt.jl:177
top-level scope

The following partial combinations work

f( (10±1)u"m" ) # works: Unitful + Measurements
quadgk(f, 0u"m", 10u"m") # works: QuadGK + Unitful
quadgk(f, 0, 10±1) # works: QuadGK + Measurements
giordano commented

Sorry, I got only now the time to look into this.

I suspect the culprit is Unitful.jl here. In this line of QuadGK, when T = typeof((10±1)u"m") we have

julia> T = typeof((10±1)u"m")
Quantity{Measurement{Float64}, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}

julia> typeof(float(real(one(T))))
Measurement{Float64}

so the Quantity information of the type is stripped away, leading to an inconsistency.

For comparison, when the number has only errors, but not units:

julia> T = typeof(10±1)
Measurement{Float64}

julia> typeof(float(real(one(T))))
Measurement{Float64}

the information of the type is preserved.

Honestly, I have no idea how to fix this in Unitful.jl (if it is indeed an issue there), but I think I have a workaround that can be applied here, I only need to make sure it doesn't break too much stuff.

DanDeepPhase commented

I wrote a bandaid fix that i wrote based on your quadgk code. It's a little brute force, but seems to work:

function quadunitmeas_result(integ::Quantity{<:Measurement},deriv,a::Measurement)
    u=unit(integ)
    integstrip=ustrip(integ)
    derivstrip=ustrip.(deriv)
    derivtuple=(1,Measurement.value.(derivstrip)...)
    meastuple=(integstrip,ustrip.(a)...)
    u*Measurements.quadgk_result(integstrip.val,derivtuple,meastuple)
end
function quadunitmeas_result(integ::Quantity,deriv,a)
    u=unit(integ)
    u*Measurements.quadgk_result(ustrip(upreferred(integ)),ustrip.(upreferred.(deriv)), ustrip.(upreferred.(a)))
    #u*Measurements.quadgk_result(ustrip((integ,deriv,a)))
end
function QuadGK.quadgk(f,a::Quantity{<:Measurement},b::Quantity;kws...)
    av=Measurements.value(a)
    integ=quadgk(f,av,b;kws...)
    deriv=-f(av)
    meas = quadunitmeas_result(integ[1],deriv,a)
    (meas,integ[2])
end
    function QuadGK.quadgk(f,a::Quantity,b::Quantity{<:Measurement};kws...)
    bv=Measurements.value(b)
    integ=quadgk(f,a,bv;kws...)
    deriv=f(bv)
    meas = quadunitmeas_result(integ[1],deriv,b)
    (meas,integ[2])
end
    function QuadGK.quadgk(f, a::Quantity{<:Measurement}, b::Quantity{<:Measurement}; kws...)
    av=Measurements.value(a)
    bv=Measurements.value(b)
    integ=quadgk(f,av,bv;kws...)
    deriv=(-f(av),f(bv))
    meas = quadunitmeas_result(integ[1],deriv,(a,b))
    (meas,integ[2])
end
giordano commented

Cool, but I think #89 is much simpler, it's just one line 😅

  • 👍1