Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Nested Types issue #75

Closed
DanBoland opened this issue Sep 8, 2020 · 3 comments · Fixed by #89
Closed

Nested Types issue #75

DanBoland opened this issue Sep 8, 2020 · 3 comments · Fixed by #89

Comments

@DanBoland
Copy link

DanBoland commented Sep 8, 2020

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 [email protected]:257[inlined]
[email protected]:257[inlined]
[email protected]: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
[email protected]:113[inlined]
#quadgk#[email protected]: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
Copy link
Member

giordano commented Feb 13, 2021

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
Copy link

DanDeepPhase commented Feb 13, 2021

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
Copy link
Member

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants