Quantity currently incompatible with QuadGK

[mikeingold]

I've been testing out DynamicQuantities as a replacement for Unitful in some of my packages and I discovered that it currently doesn't work with QuadGK's integral solver like Unitful does. MWE (issue also opened here).

julia> using DynamicQuantities  # [06fc5a27] DynamicQuantities v0.7.0

julia> using QuadGK  # [1fd47b50] QuadGK v2.8.2

julia> quadgk(t -> 5u"m/s"*t, 0u"s", 10u"s")
ERROR: MethodError: no method matching cachedrule(::Type{Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}}, ::Int64)

Closest candidates are:
  cachedrule(::Union{Type{Complex{BigFloat}}, Type{BigFloat}}, ::Integer)
   @ QuadGK C:\Users\*\.julia\packages\QuadGK\XqIlh\src\gausskronrod.jl:253
  cachedrule(::Type{T}, ::Integer) where T<:Number
   @ QuadGK C:\Users\*\.julia\packages\QuadGK\XqIlh\src\gausskronrod.jl:264

It looks like the issue probably spawns from here where Base.one(::Quantity) produces a new_quantity.

julia> typeof(one(typeof(1u"s")))
Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}

QuadGK would expect Base.one for the above example to return simply a Float64 instead. This seems consistent with the docstring for Base.one, excerpted:

If possible, one(x) returns a value of the same type as x, and one(T) returns a value of type T. However,
this may not be the case for types representing dimensionful quantities (e.g. time in days), since the
multiplicative identity must be dimensionless. In that case, one(x) should return an identity value of
the same precision (and shape, for matrices) as x.

If you want a quantity that is of the same type as x, or of type T, even if x is dimensionful, use
oneunit instead.

...

  Examples
  ≡≡≡≡≡≡≡≡≡≡

...

  julia> import Dates; one(Dates.Day(1))
  1

It also seems consistent with how Unitful handles this

julia> typeof(one(typeof(1.0u"s")))  # Unitful @u_str
Float64

Was the current definition of Base.one here intended to return the dimensioned type?

[mikeingold]

I wouldn’t mind taking a crack at a PR for this if you’re on board with changing this implementation, but I’d need to spend a little time making sure I understand the surrounding @eval for loop.

[stevengj]

Was the current definition of Base.one here intended to return the dimensioned type?

It looks like it returns a quantity with empty units, which is a correct multiplicative identity.

But returning the underlying Real type is also a correct multiplicative identity, is simpler, and allows packages like QuadGK to use one(x) to get the underlying dimensionless scalar type (for which there is currently no other cross-package API). Hence, I would tend to recommend making this change.

[MilesCranmer]

Thanks for explaining this. Yes the original implementation of one here was intentionally set to return a dimensionless quantity rather than a scalar, but I am happy to change it.

It does feel a bit strange for one(x) to not return the same type as x but if the Base docstring says so, I guess it’s fine.

[MilesCranmer]

Is that all that’s needed to get QuadGK working? We could add an integration test for it.

[MilesCranmer]

I’m happy to make the change. But, thinking out loud, I am also worried that design decisions made in Unitful (and propagated to and enforced by dependents) will influence the design of DynamicQuantities, because of network effects like this.

In particular, consider this part of the docstring:

of type T. However, this may not be the case for types representing dimensionful quantities (e.g. time in days), since the
multiplicative identity must be dimensionless.

It sounds like this only applies when the type of x stores dimension information. However, here the type of x does not store this information, and a Quantity of the exact same type can be dimensionless (and thus can be a multiplicative identity), so this does not apply. In other words, the same trick to get a dimensionless number from Unitful need not apply here.

What I’m thinking is that this usage of one was described with things like Unitful in mind, where dimensions = type. So while I’m happy to implement it to get things working in the broader ecosystem, I’m not sure it’s actually the most correct way.

But I guess we have to, if it’s a common pattern across packages.

[MilesCranmer]

@stevengj do you know other packages which assume this behavior of one? I’m just wondering if it’s a very common assumption, or if I can keep the current behavior (which feels more correct, as typeof(one(x)) == typeof(x), whereas this would break afterwards), and instead write some extensions. The practical side is that maybe you want to keep the AbstractQuantity container or AbstractDimensions type around, whereas moving to the base value type would lose that.

Maybe another option is to write a promotion rule for Quantity which attempts to convert it to its base value type, but throws a DimensionError if not dimensionless.

[MilesCranmer]

Draft in #41. It doesn't currently work, as cachedrule is only implemented for T<:Number, so it can't get to the one anyways.

Also, it doesn't call one(::Quantity), it actually calls one(::Type{Quantity}). So not sure what exactly to do for that...

It's too bad I can't do

abstract type AbstractQuantity{T,D} <: T end

which would make quantity containing a Number also a type Number.

[MilesCranmer]

Ah, so this is how Unitful.jl handles that issue:

abstract type AbstractQuantity{T,D,U} <: Number end

https://github.com/PainterQubits/Unitful.jl/blob/6fb11e7631f9b5c0897e38541ec96c0c7b3a653e/src/types.jl#L151

So I guess we could do the same, even though it feels a bit limiting.

[MilesCranmer]

Okay I fixed it. #42 is now working with QuadGK.jl

However it required the following changes. Before:

abstract type AbstractQuantity{T,D} end

one(::Quantity{T,D}) -> Quantity{T,D}
one(::Type{Quantity{T,D}}) -> Quantity{T,D}

is now

abstract type AbstractQuantity{T,D} <: Number end

one(::Quantity{T,D}) -> T
one(::Type{Quantity{T,D}}) -> T

is this actually what we want? @ChrisRackauckas or @odow I'd be eager to hear your opinion as well.

[mikeingold]

I don't have a particularly strong opinion here. The docstring for Base.one took me a while to fully parse, but I think I'm generally in favor of the interpretation where one(::Quantity{T,D}) -> T. Both versions fulfill the multiplicative identity. The docstring authors did envision that "[i]f possible, one(x) returns a value of the same type as x, and one(T) returns a value of type T" but then specify that this might not be the case for dimensionful quantities. I'm not sure if they envisioned a scenario with dimensionless quantities, though, which could still fulfill that multiplicative identity. However, the third paragraph specifies that Base.oneunit should be used when you specifically want a quantity type returned.

I kind of stumbled onto this issue after being inspired by @ChrisRackauckas over on a Julia Discourse thread. I love the premise of integrating dimensionful quantities more directly into scientific computing calculations, have had some overall good experiences with Unitful, but I think DynamicQuantities approach is more promising in the long run. I would love to have better direct support between these kinds of packages without all of the ustrip and reapplication at the interface boundaries. In any event, thanks @MilesCranmer for the fast response and the awesome package.

[stevengj]

which feels more correct, as typeof(one(x)) == typeof(x)

That's definitely not correct if the units are encoded in the type, in which case one must be a different type (even if it in the same family of types) to be a multiplicative identity. But since you don't actually do this, I guess you mean you want them to have the same supertype? But why?

Given that the type has to change, it seems a lot simpler to return the base scalar type rather than a dimensionless AbstractQuantity, though I agree that this is not required by Base.one. Moreover, there is currently no other way to get the base scalar type from dimensionful types and similar (e.g. measurements, dual numbers) AFAIK using only Base functions, so it is useful functionality to provide.