Preserve the floating-point precision of quantities in uconvert

[eliascarv]

master:

julia> using Unitful

julia> x = 100f0u"cm"
100.0f0 cm

julia> uconvert(u"m", x)
1.0f0 m

julia> x = 45f0u"°"
45.0f0°

julia> uconvert(u"rad", x)
0.7853981633974483 rad

PR:

julia> using Unitful

julia> x = 100f0u"cm"
100.0f0 cm

julia> uconvert(u"m", x)
1.0f0 m

julia> x = 45f0u"°"
45.0f0°

julia> uconvert(u"rad", x)
0.7853982f0 rad

The only drawback of this PR is that it now requires the numeric type of the quantity to implement the real function.

closes #753

[sostock]

The only drawback of this PR is that it now requires the numeric type of the quantity to implement the real function.

There are other places where real is required, e.g. the calculation of the default atol for isapprox. But I’m not sure that this means we can require it here without breaking some code.

However, even if this breaks some code, it will produce an error, so it’s not that bad (it doesn’t silently return a wrong result) and someone affected by this might open an issue.

To be on the safe side, we could add a hasmethod check (and turn the function into a @generated function so that hasmethod is only called once at compile time), like the suggested change below. But I’m not sure it is necessary. Do you know a number type that doesn’t implement real? I don’t.

[eliascarv]

Do you know a number type that doesn’t implement real? I don’t.

A Number type that is defined by another package, maybe.

This problem only occurs because the Number API is not very well defined by Julia Base.

We can assume that the Number types supported by Quantity must implement real(N) -> R<:Real and float(R) -> F<:AbstractFloat, where N is a custom Number type.

However, the solution you recommended seems good.
What do you think, @sostock?

[sostock]

We can assume that the Number types supported by Quantity must implement real(N) -> R<:Real and float(R) -> F<:AbstractFloat, where N is a custom Number type.

This is unfortunately not the case for the types from CliffordNumbers.jl (which supports Unitful numbers via a package extension). They implement real and float, but return Clifford numbers again, which are not <: Real:

julia> k = KVector{1, VGA(3)}(1, 2, 3)
3-element KVector{1, VGA(3), Int64}:
1e₁ + 2e₂ + 3e₃

julia> real(k)
3-element KVector{1, VGA(3), Int64}:
1e₁ + 2e₂ + 3e₃

julia> real(k) isa Real
false

julia> float(k)
3-element KVector{1, VGA(3), Float64}:
1.0e₁ + 2.0e₂ + 3.0e₃

julia> float(k) isa AbstractFloat
false

This means that CliffordNumbers would not benefit from this change (a Float32-based Clifford number would still be widened to a Float64-based Clifford number on certain unit conversions).

To make this non-widening behavior easy to opt-in, we could add a new function that packages can overload for their types. The fallback would be this:

function floattype(::Type{T}) where T
    try # Use try-catch instead of hasmethod because real(T) fallback method exists but might throw an error
        float(real(T))
    catch
        Nothing
    end
end

Packages like CliffordNumbers.jl can then add their own Unitful.floattype methods. In convfact, we check whether floattype(T) <: AbstractFloat and, if that is the case, use its value.

What do you think?

[eliascarv]

@sostock the Unitful.floattype seems like a good idea.
But maybe floattype could be written like this:

floattype(::Type{T}) where {T<:AbstractFloat} = T
floattype(::Type{T}) where {T<:Complex} = floattype(real(T))
floattype(::Type{T}) where {T<:Number} = Float64 # or nothing, but Float64 seems like a good default

What do you think?

[sostock]

But maybe floattype could be written like this:

floattype(::Type{T}) where {T<:AbstractFloat} = T
floattype(::Type{T}) where {T<:Complex} = floattype(real(T))
floattype(::Type{T}) where {T<:Number} = Float64 # or nothing, but Float64 seems like a good default

That would mean floattype(Quaternion{Float32}) == Float64, so Quaternions.jl would have to implement their own floattype method even though float(real(T)) would work correctly for Quaternions.

[eliascarv]

@sostock you're absolutely right. I think the most viable solution is actually the try catch. Maybe:

function floattype(::Type{T}) where {T<:Number}
    try
        F = float(real(T))
        F isa AbstractFloat ? F : Float64 # or nothing
    catch
        Float64 # or nothing
    end
end

Since this function will run at compile time, I think the cost of the try catch will not be a problem.

[cstjean]

I think the cost of the try catch will not be a problem.

Famous last words!

I would at least benchmark and compare uconvert.(u"...", rand(100000) .* u"...")) using

• Current master
• try/catch
• if applicable(...)

Also, is this package tested at all against the autodiff packages? I don't know how fragile those are, but try/catch feels risky.

[sostock]

I have already benchmarked some cases and will benchmark some more when we have decided which version to implement. I have actually found a slight performance regression (that is fortunately easy to fix), in the original version of the PR (the one without the try … catch):

julia> using BenchmarkTools, 

julia> @btime uconvert(u"rad", $(big(45.0)u"°"));
  109.262 ns (2 allocations: 104 bytes) # v1.21.1: conversion factor is a Float64
  122.539 ns (4 allocations: 208 bytes) # this PR: conversion factor is converted to BigFloat

It is unfortunate that the conversion factor in this case only has Float64 precision, but that is another issue.

Regarding autodiff packages, I have to admit that I don’t have any experience with any of them, and they are not tested as part of our test suite. Since try … catch will only be used in a @generated function (and not in the expression it returns), and if it fails we fall back to the old behavior, I don’t think it would break anything, but it’s better to test it. Easiest would be to look whether the autodiff packages have tests for Unitful quantities and run their test suite against this PR.

However, I also thought about another way of implementing this that shouldn’t require other packages to implement floattype to opt into this behavior: Whenever convfact would return a floating-point number, it could instead return an AbstractIrrational. We would create a new type UnitConversionFactor <: AbstractIrrational that wraps the resulting conversion factor and implements all the required methods. Of course, that could cause problems for autodiff (and maybe other packages) as well, it will need to be tested extensively. But I think it would be the more elegant way, and could be a first step towards arbitrary-precision conversion factors.

[eliascarv]

Great idea, @sostock. I think I can implement this idea in this same PR.
But I have a question: Will UnitConversionFactor be used for integer and rational conversion factors? Or only for float conversion factors?
In other words: convfactor(s, t) -> UnitConversionFactor or convfactor(s, t) -> UnitConversionFactor | Rational{Int} | Int
Maybe the second option is the correct answer?

[sostock]

Will UnitConversionFactor be used for integer and rational conversion factors? Or only for float conversion factors?
Maybe the second option is the correct answer?

Yes, I would only use it for float conversion factors. The main (or only) purpose of the AbstractIrrational type is picking an appropriate floating-point precision, so using it also for Integer and Rational will make things more complicated (I think).

[sostock]

Suggested change

	Base.hash(x::UnitConversionFactor, h::UInt) = 3 * objectid(x) - h
	Base.hash(::UnitConversionFactor{x}, h::UInt) where {x} = hash(x, h)

Since UnitConversionFactor is (for now) always exactly equal to a Float64, it should hash the same. We should also implement ==(::UnitConversionFactor, ::Real) and ==(::Real, ::UnitConversionFactor). However, since this type is only used internally, it probably doesn’t matter too much.

[eliascarv]

@sostock adjusted. Do you think the == with Real is necessary? If so, I can add the methods and tests.

[sostock]

It probably isn’t necessary, since == or hash is likely never used on this type. But I still think hash and isequal (which falls back to ==) should be consistent. So, I would either change back hash to the old implementation (where the hashes are never equal to those of a Float64) or add the == methods.

I think we could actually remove the hash and == methods and they would still be consistent using the fallback methods.

[eliascarv]

@sostock ready for a new review.

[juliohm]

@sostock appreciate if you could take a look into this.