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.

[MilesCranmer]

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?

@stevengj I'm not sure if you've read the README so apologies if I'm overexplaining here: the underlying motivation for DynamicQuantities.jl is precisely this. DynamicQuantities stores units in the value rather than the type. This avoids type instabilities and overspecialized methods. Even in type stable situations, there is not much of a slowdown either – see example on README. Chris explains the benefits better than I can here.

Thus, in DynamicQuantities, we can both have one(q) return a multiplicative identity for q, AND have that return value be same type as the input. There is not really any reason for one(::Quantity) to return a value of a different type, other than if we are simply trying to maintain compatibility with Unitful.jl (which is very fair!). Does this help explain things?

This is my main issue here. We no longer need to do this, because with DynamicQuantities, you can have the same exact type of quantity be a multiplicative identity.

But in any case I am happy to implement it. I think that Unitful.jl compatibility, where possible, should be prioritized.

[gaurav-arya]

Cross-referencing https://github.com/SymbolicML/DynamicQuantities.jl/pull/42#issuecomment-1726343281: changing one has the potential to make a function like Base.prod type unstable. This wasn't a problem with Unitful because something like a multiplicative reduction on an array was inherently type unstable anyway. However, DynamicQuantities tries very hard not to introduce type instabilities that the corresponding unitless version does not have. Thus, it seems like an alternative approach to one(T) (via an extension or otherwise?) might be more appropriate for interoperability with QuadGK.jl?

[odow]

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

Where possible, I'd keep the one(::Quantity{T,D}) -> Quantity{T,D}. There are so many places in the standard library that have implicit assumptions that one(T)::T, but most are hidden behind promotion rules so that they mostly work.

(Although JuMP's main problem is that VariableRef does not have a type-stable multiplicative or additive identity, because x::VariableRef + zero(::VariableRef)::AffExpr is an AffExpr not a VariableRef. Since the DynamicQuantities case is the other way around, either is probably fine.)

[MilesCranmer]

Thanks for bringing up these points. It's very useful to know that one(T)::T is an assumption in the standard library, and about the potential type instabilities. I think those points rule out changing the behavior of one(::Type{Quantity{T,D}}), so we can leave it as-is.

With this in mind maybe we could:

1. Write an extension for QuadGK.jl (either in DynamicQuantities or QuadGK, whichever makes the most sense)
2. Encourage packages to use a function like:

"Get the first parametric type, if it exists"
function base_type(::Type{T}) where {T}
    params = T isa UnionAll ? T.body.parameters : T.parameters
    return length(params) == 0 ? T : params[1]
end

to generically unwrap types like quantities and dual numbers, rather than using one, which seems like incorrect behavior?

---

Aside: it would be nice if something like this function was in the standard library so packages could overload its behavior in case the base type is not listed first. Rather than needing to rely on ambiguous behavior from Base.one.

Default behavior:

Quantity{Float64} -> Float64
Array{Float16,2} -> Float16
Float32 -> Float32
ComplexF64 -> Float64
Array{ComplexF64,1} -> ComplexF64
Dual{Int64} -> Int64

with failure cases like:

julia> base_type(typeof(StaticArrays.SArray{Tuple{32,3}}(randn(32, 3))))
Tuple{32, 3}

that would need to be overloaded... (or with a custom method for AbstractArray)

[odow]

The main things that would be good to have tests for are various combinations of SparseArrays and the LinearAlgebra matrix types. it mainly comes up when the matrix has a structural assumption that some elements are zero or one.

[stevengj]

There are so many places in the standard library that have implicit assumptions that one(T)::T

Really? If that's true, a lot of them are probably bugs, since that's only supposed to be guaranteed for oneunit(T)::T. This is the whole point of why separate one and oneunit functions exist.

(Most of the places where Base assumes one(T)::T that I'm aware of are for narrower types like T<:Integer, which can't really be dimensionful.)

That being said, if you store the units in the value rather than the type (to make things like prod type-stable), I can see why you'd prefer to return the same type from one. (prod does not assume that one(T)::T ... It's just type-unstable if that's not the case, which is unavoidable for things like Unitful that store the units in the type.)

[stevengj]

Aside: it would be nice if something like this function was in the standard library

I agree.

In the short term I suppose it could be provided by a little package that defines a function base_numeric_type (or whatever), which is depended on (and overloaded) by all of the unitful types, Measurements, ForwardDiff, etcetera.

[MilesCranmer]

I can make a package for this. Just to confirm, would you need Dual{ComplexF64} to return Float64 or ComplexF64?

Also, ComplexF64 -> ComplexF64 or Float64?

Also would you recommend making ext for special cases, or just ask the specific package to make an ext?

[stevengj]

Complex64 is fine, since we already have real(T) to get the real type from a complex (or quaternion, etc) numeric type.

Not sure what you are asking about with ext. I don't think this should be a package extension — it should just be a tiny package that things like QuadGK, DynamicQuantities, Measurements, etcetera can depend on unconditionally.)

[MilesCranmer]

Got it, thanks.

The ext question was more: say I want to accelerate the transition to use this theoretical package. For special cases where the above function does not work for a particular wrapper (?), I was thinking whether it makes more sense to define those special cases in this package within ext, or let things like Measurements.jl write their own interface. (Not specifically Measurements.jl, as it would already work based on the generic function above — just whatever libraries need special behavior)

[MilesCranmer]

Actually I can't think of any special cases after adding a mapping AbstractArray{T} -> T. So it's probably fine.

I created https://github.com/SymbolicML/BaseType.jl and gave you maintain access @stevengj. So the idea is that the base_numeric_type in that package can handle all cases where one was used in the past for inferring the base numeric type, as well as new edge cases as discussed in this thread. And once that function enters regular usage we can lobby for moving it to stdlib.

[stevengj]

Some points:

1. This package should not depend on Measurements, Unitful, etcetera. It should depend on nothing. Instead, the dependencies should go the other way around — Measurements etcetera should add a dependency on BaseType.jl, and extend base_numeric_type as needed. That way, a package like QuadGK can depend on BaseType.jl without pulling in zillions of other packages.

2. I don't think it should work work on containers (Array, Set, etcetera), only for ::Number. We already have eltype for containers, and it's conceptually a different sort of function.

3. Probably the package should have a more specific name, like BaseNumericType, and probably should go into JuliaMath rather than SymbolicML.

[MilesCranmer]

Some points:

1. This package should not depend on Measurements, Unitful, etcetera. It should depend on nothing. Instead, the dependencies should go the other way around — Measurements etcetera should add a dependency on BaseType.jl, and extend base_numeric_type as needed. That way, a package like QuadGK can depend on BaseType.jl without pulling in zillions of other packages.

Don't worry, this wasn't what I was suggesting. If I would do something in this direction, each interface would go in ext and thus not add additional dependencies. But it's obviously not a great solution so I'll avoid it.

2. I don't think it should work work on containers (Array, Set, etcetera), only for ::Number. We already have eltype for containers, and it's conceptually a different sort of function.

Good point, I'll remove them from the examples.

3. Probably the package should have a more specific name, like BaseNumericType,

Not sure about this one:

• Con: I could see myself adding other functions that get base types for different kinds of wrappers – including non-numerical ones.
• Pro: But, I wouldn't want to have this tiny interface package try to do too much.
• Con: I don't want to maintain several separate packages for different type extraction utilities. One package would be easier.
• Pro: package is the same name as the method, easy to remember

Also there's no other Base*Type so there wouldn't be ambiguity.

That being said I don't feel strongly about this so if you would prefer BaseNumericType.jl, I will go ahead and change it.

and probably should go into JuliaMath rather than SymbolicML.

If someone can move it there and give me write or maintain access, then that sounds good to me!

[MilesCranmer]

FYI the countdown for the package to be registered ends tomorrow, so if you would prefer BaseNumericType.jl, let me know soon.

[mikeingold]

Update: looks like this Issue is still valid but now errors differently. Integration of functions that output DynamicQuantities dimensionful types works on a primitive/float-valued domain but not on a dimensionful domain.

# Julia v1.10.2, Windows x86-64
# Temp environment with only:
#   DynamicQuantities v0.12.3
#   QuadGK v2.9.4

julia> using DynamicQuantities, QuadGK

# Integrate a function with dimensionful output over a scalar-valued domain
julia> quadgk(t -> 5u"m/s"*t, 0, 10)
(250.0 m s⁻¹, 0.0 m s⁻¹)

# Integrate a function with dimensionful output over a dimensionful domain
julia> quadgk(t -> 5u"m/s"*t, 0u"s", 10u"s")
MethodError: no method matching kronrod(::Type{Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}}, ::Int64)

Closest candidates are:
  kronrod(::Any, ::Integer, ::Real, ::Real; rtol, quad)
   @ QuadGK C:\Users\mikei\.julia\packages\QuadGK\OtnWt\src\weightedgauss.jl:90
  kronrod(::Type{T}, ::Integer) where T<:AbstractFloat
   @ QuadGK C:\Users\mikei\.julia\packages\QuadGK\OtnWt\src\gausskronrod.jl:316
  kronrod(::AbstractMatrix{<:Real}, ::Integer, ::Real, ::Pair{<:Tuple{Real, Real}, <:Tuple{Real, Real}})
   @ QuadGK C:\Users\mikei\.julia\packages\QuadGK\OtnWt\src\gausskronrod.jl:390
  ...

Stacktrace:
 [1] macro expansion
   @ C:\Users\mikei\.julia\packages\QuadGK\OtnWt\src\gausskronrod.jl:564 [inlined]
 [2] _cachedrule
   @ C:\Users\mikei\.julia\packages\QuadGK\OtnWt\src\gausskronrod.jl:564 [inlined]
 [3] cachedrule
   @ C:\Users\mikei\.julia\packages\QuadGK\OtnWt\src\gausskronrod.jl:569 [inlined]
 [4] do_quadgk(f::var"#1#2", s::Tuple{Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}, Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64, nrm::typeof(LinearAlgebra.norm), segbuf::Nothing)
   @ QuadGK C:\Users\mikei\.julia\packages\QuadGK\OtnWt\src\adapt.jl:7
 [5] (::QuadGK.var"#50#51"{Nothing, Nothing, Int64, Int64, typeof(LinearAlgebra.norm), Nothing})(f::Function, s::Tuple{Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}, Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}}, ::Function)
   @ QuadGK C:\Users\mikei\.julia\packages\QuadGK\OtnWt\src\adapt.jl:253
 [6] handle_infinities(workfunc::QuadGK.var"#50#51"{Nothing, Nothing, Int64, Int64, typeof(LinearAlgebra.norm), Nothing}, f::var"#1#2", s::Tuple{Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}, Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}})
   @ QuadGK C:\Users\mikei\.julia\packages\QuadGK\OtnWt\src\adapt.jl:145
 [7] quadgk(::Function, ::Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}, ::Vararg{Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}}; atol::Nothing, rtol::Nothing, maxevals::Int64, order::Int64, norm::Function, segbuf::Nothing)
   @ QuadGK C:\Users\mikei\.julia\packages\QuadGK\OtnWt\src\adapt.jl:252
 [8] quadgk(::Function, ::Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}, ::Vararg{Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}})
   @ QuadGK C:\Users\mikei\.julia\packages\QuadGK\OtnWt\src\adapt.jl:250
 [9] top-level scope
   @ REPL[5]:1

There's a kronrod method for <:AbstractFloat types here. I don't see any other two-argument methods for other non-float types. Do Unitful.jl Quantity's sub-type from AbstractFloat? I don't see how they would even work otherwise.

[MilesCranmer]

In principle we can add a FloatQuantity (to complement the existing Quantity, RealQuantity, GenericQuantity), but it might increase loading times by a factor.

[MilesCranmer]

Oh but right, Unitful is also <: Number.. hmm..