Base.isapprox fails for Vector{Quantity}

[mikeingold]

Extending this from a comment in #76, I reduced the issue to a MWE. The issue seems to occur specifically when running Base.isapprox on a Vector{Quantity}.

using DynamicQuantities

a = 1.0u"m"
b = 2.0u"m"

a ≈ b
# Works ok: returns False

[a] ≈ [b]
#=
ERROR: Cannot create an additive identity for a `UnionAbstractQuantity` type, as the dimensions are unknown. Please use `zero(::UnionAbstractQuantity)` instead.
Stacktrace:
 [1] error(s::String)
   @ Base .\error.jl:35
 [2] zero(::Type{Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}})
   @ DynamicQuantities C:\Users\mikei\.julia\packages\DynamicQuantities\jkfTz\src\utils.jl:280
 [3] real(T::Type)
   @ Base .\complex.jl:120
 [4] rtoldefault(x::Type{Quantity{Float64, Dimensions{…}}}, y::Type{Quantity{Float64, Dimensions{…}}}, atol::Int64)
   @ Base .\floatfuncs.jl:347
 [5] isapprox(x::Vector{Quantity{Float64, Dimensions{…}}}, y::Vector{Quantity{Float64, Dimensions{…}}})
   @ LinearAlgebra C:\Users\mikei\.julia\juliaup\julia-1.10.0+0.x64.w64.mingw32\share\julia\stdlib\v1.10\LinearAlgebra\src\generic.jl:1785
 [6] top-level scope
   @ REPL[8]:1
=#

[mikeingold]

Looks like I mis-stated in #76: the method isapprox(::AbstractArray, ::AbstractArray) is implemented under LinearAlgebra.jl here, currently:

function isapprox(x::AbstractArray, y::AbstractArray;
    atol::Real=0,
    rtol::Real=Base.rtoldefault(promote_leaf_eltypes(x),promote_leaf_eltypes(y),atol),
    nans::Bool=false, norm::Function=norm)
    d = norm(x - y)
    if isfinite(d)
        return d <= max(atol, rtol*max(norm(x), norm(y)))
    else
        # Fall back to a component-wise approximate comparison
        return all(ab -> isapprox(ab[1], ab[2]; rtol=rtol, atol=atol, nans=nans), zip(x, y))
    end
end

For the MWE, the rtol arg currently resolves to

Base.rtoldefault(
    Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}},
    Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}},
    atol
)

[MilesCranmer]

Hm, does it work if you do RealQuantity instead of Quantity? That would let you pass it to Real arguments. RealQuantity identical to Quantity except it is <:Real rather than <:Number – the purpose is for passing values to such functions.

e.g.,

RealQuantity(1.0u"m")

[MilesCranmer]

I also realised our current implementation of isapprox doesn't actually do anything special for atol – it just assumes the bound is going to be computed on the stripped quantity:

DynamicQuantities.jl/src/utils.jl

Lines 224 to 236 in 9681769

	function Base.isapprox(l::$type, r::$type; kws...)
	l, r = promote_except_value(l, r)
	dimension(l) == dimension(r) || throw(DimensionError(l, r))
	return isapprox(ustrip(l), ustrip(r); kws...)
	end
	function Base.isapprox(l::$base_type, r::$type; kws...)
	iszero(dimension(r)) || throw(DimensionError(l, r))
	return isapprox(l, ustrip(r); kws...)
	end
	function Base.isapprox(l::$type, r::$base_type; kws...)
	iszero(dimension(l)) || throw(DimensionError(l, r))
	return isapprox(ustrip(l), r; kws...)
	end

So maybe we want to have custom logic for atol?

What is the assumption in Unitful? We can just match the existing API on this sort of thing.

[mikeingold]

Hm, does it work if you do RealQuantity instead of Quantity? That would let you pass it to Real arguments. RealQuantity identical to Quantity except it is <:Real rather than <:Number – the purpose is for passing values to such functions.

This is definitely an option, but it’s also pretty natural in practice to write something like [1.0u”m/s”, 2.0u”m/s”] and anticipate that Base.isapprox will understand how to compare it.

[mikeingold]

#What is the assumption in Unitful? We can just match the existing API on this sort of thing.

It looks like they defined a bunch of Base.isapprox methods here:

https://github.com/PainterQubits/Unitful.jl/blob/a09bdaac255d456e6315b77979a48956aca1d3a3/src/quantities.jl#L248-L295

A call to @which [a] ≈ [b] under Unitful indicates that it calls this one, which looks like a modified copy of the one from LinearAlgebra above.

[mikeingold]

It seems like the most robust way to handle isapprox with tolerances on DQ types would be to just define specific methods for them. I'm taking a quick pass now to see what we would need to define.

For defining default rtol's, the core of it seems to come down to the following, where floats get set to sqrt(eps) and everything else gets set to zero.

rtoldefault(::Type{T}) where {T<:AbstractFloat} = sqrt(eps(T))
rtoldefault(::Type{<:Real}) = 0

These should be reasonably simple to reimplement to something like:

rtoldefault(x::UnionAbstractQuantity{T,D}) where {T<:AbstractFloat,D} = sqrt(eps(T)) * x.dimensions
rtoldefault(x::UnionAbstractQuantity{T,D}) where {T<:Real,D} = zero(x)

[MilesCranmer]

That sounds good to me!! Or even something like

rtoldefault(q::UnionAbstractQuantity) = new_quantity(typeof(q), rtoldefault(ustrip(q)), dimension(q))

to be even more generic, as it would preserve the quantity type and also take the underlying rtoldefault for T.

[MilesCranmer]

In this case I think we can simply add :rtoldefault to this list:

DynamicQuantities.jl/src/math.jl

Lines 193 to 194 in 9681769

	:float, :abs, :real, :imag, :conj, :adjoint, :unsigned,
	:nextfloat, :prevfloat, :transpose, :significand

since it has the same structure

[mikeingold]

Sounds good to me. There’s definitely going to be a little more involved since the default isapprox enforces tol::Real for both in the args list, but it doesn’t seem too complicated to re-define overall. I haven’t had time to dig into it again since last night, but can probably take another crack at it later.

[mikeingold]

Alright, I started working more on an implementation for this and realized I needed to do some disambiguation in my head about what exactly isapprox actually answers. I'm putting together a PR to get the tolerance kwargs working in general and for proper isapprox on vectors.

For the case of isapprox(a, b) where a and b are both numbers, with equivalent quantity-dimensions, the result should be a Boolean true/false according to whether some distance measure (norm) between a and b is within some absolute tolerance (atol) or relative tolerance (rtol) threshold.

For the case of isapprox(u, v) where u and v are both $n$-dimensional vectors, the result should be whether or not the norm in $n$-space between these two vectors is within the specified tolerances. The element-wise pairs between u and v must have equivalent quantity-dimensions in order for norm to be well-defined. This operation isapprox(u, v) is not simply an element-wise comparison. That would be achieved by the user instead using something like all(isapprox.(u, v)).

[mikeingold]

This seems pretty simple when written out, but what tripped me up was trying to think about something like an isapprox(u, v) where u and v are state-space vectors. For example, let u = [1.0u"m", 2.0u"m/s"] and v = [2.0u"m", 3.0u"m/s"]. I don't have a strong enough maths background to say for certain whether a single-valued distance measure truly exists in such a state space (in units $\sqrt{m^2+(m/s)^2}$?), but it doesn't seem directly representable in Unitful or DynamicQuantities.

I think the most practical answer in this kind of situation is to require that the quantity-dimensions of the vector elements be equivalent in order to compute isapprox(u, v). In a case where the element-wise pairs (e.g. zip(u,v)) have the same quantity-dimensions but not all of these units are equivalent along the vectors, we might want to detect this condition and throw an error to help explain the conundrum (like "no well-defined/calculable norm available, consider direct element-wise comparison"). In a case where element-wise pairs don't even have the same units, though, I'm undecided but lean toward throwing an error (like "can't compute norm due to dimensional incompatibility"); returning false implies that norm(u-v) >= tol but this norm isn't even calculable.

[MilesCranmer]

In a case where element-wise pairs don't even have the same units, though, I'm undecided but lean toward throwing an error (like "can't compute norm due to dimensional incompatibility")

I think this is a good idea. You could do something like

(allequal(dimension.(u)) && allequal(dimension.(v)) && dimension(first(u)) == dimension(first(v))) ||
    throw(DimensionError(u, v))

Basically whatever prevents the user from shooting themselves in the foot here!

Actually if you just compute norm of a vector that isn't dimensionally consistent, it will already throw an error without this change.

Even if that means basically throwing an error on isapprox(u, v) and forcing the user to call ustrip themselves, I think that's an option too.

[mikeingold]

Update: I’m still working on this as time permits.

I think I’ve gotten the dimensionful atol implemented; it required adding a helper function to inspect kwargs, compare the dimension of an optional atol and return a ustrip’d version of it if found good. This is probably going to be considered a breaking change and fails a number of current unit tests in the Unitful extension since those use isapprox on quantities with dimensionless atols.

Upon further reflection, I don’t think rtol is supposed to be dimensionful. Relative to me implies a ratio, which should be unitless, so maybe ok as-is?

I’ve got a rough sketch of the vector comparison version written that I suspect will work fine once I’ve gotten the regular version nailed down and any lessons learned incorporated.

[mikeingold]

Fixed with merge of PR #114