add Test.GenericDimensionful number type

NEWS.md

@@ -39,11 +39,16 @@ New library features
 --------------------
 
 * The optional keyword argument `context` of `sprint` can now be set to a tuple of `:key => value` pairs to specify multiple attributes. ([#39381])
+* `atol` argument of `isapprox` now accepts arbitrary `Number` types (rather than `Real`)
+  in order to support dimensionful values (though a nonzero `atol` must still support `max`
+  like a real number) ([#39852]).
 
 Standard library changes
 ------------------------
 
 * `count` and `findall` now accept an `AbstractChar` argument to search for a character in a string ([#38675]).
+* `Test` now exports a `GenericDimensionful` number type that can be used to test whether
+  code works with dimensionful types (like those in the Unitful.jl package) ([#39852]).
 * `range` now supports the `range(start, stop)` and `range(start, stop, length)` methods ([#39228]).
 * `range` now supports `start` as an optional keyword argument ([#38041]).
 * `islowercase` and `isuppercase` are now compliant with the Unicode lower/uppercase categories ([#38574]).

base/floatfuncs.jl

@@ -223,16 +223,16 @@ end
 
 # isapprox: approximate equality of numbers
 """
-    isapprox(x, y; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps, nans::Bool=false[, norm::Function])
+    isapprox(x, y; atol::Number=0, rtol::Real=atol>0 ? 0 : √eps, nans::Bool=false[, norm::Function])
 
 Inexact equality comparison: `true` if `norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))`. The
 default `atol` is zero and the default `rtol` depends on the types of `x` and `y`. The keyword
 argument `nans` determines whether or not NaN values are considered equal (defaults to false).
 
-For real or complex floating-point values, if an `atol > 0` is not specified, `rtol` defaults to
+For real or complex floating-point values, if an `atol ≠ 0` is not specified, `rtol` defaults to
 the square root of [`eps`](@ref) of the type of `x` or `y`, whichever is bigger (least precise).
 This corresponds to requiring equality of about half of the significand digits. Otherwise,
-e.g. for integer arguments or if an `atol > 0` is supplied, `rtol` defaults to zero.
+e.g. for integer arguments or if an `atol ≠ 0` is supplied, `rtol` defaults to zero.
 
 The `norm` keyword defaults to `abs` for numeric `(x,y)` and to `LinearAlgebra.norm` for
 arrays (where an alternative `norm` choice is sometimes useful).
@@ -257,6 +257,10 @@ but an absurdly large tolerance if `x` is the
     Passing the `norm` keyword argument when comparing numeric (non-array) arguments
     requires Julia 1.6 or later.
 
+!!! compat "Julia 1.7"
+    Passing a `Number` value for `atol` rather than `Real` (which is useful for
+    dimensionful number types) requires Julia 1.7 or later.
+
 # Examples
 ```jldoctest
 julia> 0.1 ≈ (0.1 - 1e-10)
@@ -276,11 +280,15 @@ true
 ```
 """
 function isapprox(x::Number, y::Number;
-                  atol::Real=0, rtol::Real=rtoldefault(x,y,atol),
+                  atol::Number=0, rtol::Real=rtoldefault(x,y,atol),
                   nans::Bool=false, norm::Function=abs)
-    x == y || (isfinite(x) && isfinite(y) && norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))) || (nans && isnan(x) && isnan(y))
+    x == y || (isfinite(x) && isfinite(y) && _isapprox_small(norm(x-y), atol, rtol*max(norm(x), norm(y)))) || (nans && isnan(x) && isnan(y))
 end
 
+# check if d is small compared to atol and rtolx (= rtol*scale),
+# ignoring the mismatch in units (if any) between d and atol only when atol is zero.
+_isapprox_small(d, atol, rtolx) = d <= (iszero(atol) ? rtolx : max(atol, rtolx))
+
 """
     isapprox(x; kwargs...) / ≈(x; kwargs...)
 
@@ -301,9 +309,10 @@ This is equivalent to `!isapprox(x,y)` (see [`isapprox`](@ref)).
 # default tolerance arguments
 rtoldefault(::Type{T}) where {T<:AbstractFloat} = sqrt(eps(T))
 rtoldefault(::Type{<:Real}) = 0
-function rtoldefault(x::Union{T,Type{T}}, y::Union{S,Type{S}}, atol::Real) where {T<:Number,S<:Number}
+rtoldefault(::Type{T}) where {T<:Number} = rtoldefault(typeof(real(one(T)))) # strip dimensions if any
+function rtoldefault(x::Union{T,Type{T}}, y::Union{S,Type{S}}, atol::Number) where {T<:Number,S<:Number}
     rtol = max(rtoldefault(real(T)),rtoldefault(real(S)))
-    return atol > 0 ? zero(rtol) : rtol
+    return iszero(atol) ? rtol : zero(rtol)
 end
 
 # fused multiply-add

base/math.jl

@@ -611,7 +611,7 @@ function _hypot(x, y)
 
     # Return Inf if either or both inputs is Inf (Compliance with IEEE754)
     if isinf(ax) || isinf(ay)
-        return oftype(axu, Inf)
+        return typeof(axu)(Inf)
     end
 
     # Order the operands
@@ -661,7 +661,7 @@ _hypot(x::ComplexF16, y::ComplexF16) = Float16(_hypot(ComplexF32(x), ComplexF32(
 function _hypot(x...)
     maxabs = maximum(abs, x)
     if isnan(maxabs) && any(isinf, x)
-        return oftype(maxabs, Inf)
+        return typeof(maxabs)(Inf)
     elseif (iszero(maxabs) || isinf(maxabs))
         return maxabs
     else

stdlib/LinearAlgebra/src/generic.jl

@@ -478,7 +478,7 @@ function generic_norm1(x)
     (v, s) = iterate(x)::Tuple
     av = float(norm(v))
     T = typeof(av)
-    sum::promote_type(Float64, T) = av
+    sum = av*1.0 # promotes to at least Float64 with * 1.0
     while true
         y = iterate(x, s)
         y === nothing && break
@@ -495,11 +495,11 @@ norm_sqr(x::Union{T,Complex{T},Rational{T}}) where {T<:Integer} = abs2(float(x))
 
 function generic_norm2(x)
     maxabs = normInf(x)
-    (maxabs == 0 || isinf(maxabs)) && return maxabs
+    (iszero(maxabs) || isinf(maxabs)) && return maxabs
     (v, s) = iterate(x)::Tuple
     T = typeof(maxabs)
-    if isfinite(length(x)*maxabs*maxabs) && maxabs*maxabs != 0 # Scaling not necessary
-        sum::promote_type(Float64, T) = norm_sqr(v)
+    if isfinite(length(x)*maxabs*maxabs) && !iszero(maxabs*maxabs) # Scaling not necessary
+        sum = norm_sqr(v)*1.0 # promotes to at least Float64 with * 1.0
         while true
             y = iterate(x, s)
             y === nothing && break
@@ -525,30 +525,32 @@ function generic_normp(x, p)
     (v, s) = iterate(x)::Tuple
     if p > 1 || p < -1 # might need to rescale to avoid overflow
         maxabs = p > 1 ? normInf(x) : normMinusInf(x)
-        (maxabs == 0 || isinf(maxabs)) && return maxabs
+        (iszero(maxabs) || isinf(maxabs)) && return maxabs
         T = typeof(maxabs)
+        maxabsd = maxabs*1.0 # promotes to at least Float64 with * 1.0
     else
         T = typeof(float(norm(v)))
+        maxabsd = oneunit(T)*1.0 # promotes to at least Float64 with * 1.0
     end
-    spp::promote_type(Float64, T) = p
-    if -1 <= p <= 1 || (isfinite(length(x)*maxabs^spp) && maxabs^spp != 0) # scaling not necessary
-        sum::promote_type(Float64, T) = norm(v)^spp
+    pinv = p isa Integer ? 1//p : 1.0*inv(p)
+    if -1 <= p <= 1 || (isfinite(length(x)*maxabsd^p) && !iszero(maxabsd^p)) # scaling not necessary
+        sum = (norm(v) * 1.0)^p # promotes to at least Float64 with * 1.0
         while true
             y = iterate(x, s)
             y === nothing && break
             (v, s) = y
-            sum += norm(v)^spp
+            sum += (norm(v) * 1.0)^p
         end
-        return convert(T, sum^inv(spp))
+        return convert(T, sum^pinv)
     else # rescaling
-        sum = (norm(v)/maxabs)^spp
+        sum = (norm(v)/maxabsd)^p
         while true
             y = iterate(x, s)
             y === nothing && break
             (v, s) = y
-            sum += (norm(v)/maxabs)^spp
+            sum += (norm(v)/maxabsd)^p
         end
-        return convert(T, maxabs*sum^inv(spp))
+        return convert(T, maxabsd*sum^pinv)
     end
 end
 
@@ -1654,12 +1656,12 @@ promote_leaf_eltypes(x::Union{AbstractArray,Tuple}) = mapreduce(promote_leaf_elt
 # Supports nested arrays; e.g., for `a = [[1,2, [3,4]], 5.0, [6im, [7.0, 8.0]]]`
 # `a ≈ a` is `true`.
 function isapprox(x::AbstractArray, y::AbstractArray;
-    atol::Real=0,
+    atol::Number=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)))
+        return Base._isapprox_small(d, 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))

stdlib/LinearAlgebra/src/uniformscaling.jl

@@ -328,18 +328,18 @@ isequal(A::AbstractMatrix, J::UniformScaling) = false
 isequal(J::UniformScaling, A::AbstractMatrix) = false
 
 function isapprox(J1::UniformScaling{T}, J2::UniformScaling{S};
-            atol::Real=0, rtol::Real=Base.rtoldefault(T,S,atol), nans::Bool=false) where {T<:Number,S<:Number}
+            atol::Number=0, rtol::Real=Base.rtoldefault(T,S,atol), nans::Bool=false) where {T<:Number,S<:Number}
     isapprox(J1.λ, J2.λ, rtol=rtol, atol=atol, nans=nans)
 end
 function isapprox(J::UniformScaling, A::AbstractMatrix;
-                  atol::Real = 0,
+                  atol::Number = 0,
                   rtol::Real = Base.rtoldefault(promote_leaf_eltypes(A), eltype(J), atol),
                   nans::Bool = false, norm::Function = norm)
     n = checksquare(A)
     normJ = norm === opnorm             ? abs(J.λ) :
             norm === LinearAlgebra.norm ? abs(J.λ) * sqrt(n) :
                                           norm(Diagonal(fill(J.λ, n)))
-    return norm(A - J) <= max(atol, rtol * max(norm(A), normJ))
+    return Base._isapprox_small(norm(A - J), atol, rtol * max(norm(A), normJ))
 end
 isapprox(A::AbstractMatrix, J::UniformScaling; kwargs...) = isapprox(J, A; kwargs...)
 

stdlib/LinearAlgebra/test/generic.jl

@@ -504,4 +504,13 @@ end
     @test condskeel(A) ≈ condskeel(A, [8,8,8])
 end
 
+@testset "dimensionful arrays" begin
+    for p in (1,2,3,Inf), x in ([0,0,0], [1e300,2e300], [1,2,Inf])
+        @test GenericDimensionful(norm(x, p)) == norm(GenericDimensionful.(x), p)
+    end
+    @test dot(GenericDimensionful.([1,2,3]), [4,5,6]) == GenericDimensionful(32)
+    @test dot(GenericDimensionful.([1,2,3]), GenericDimensionful.([4,5,6])) == GenericDimensionful{2}(32)
+    @test GenericDimensionful.([1,0]) ≈ GenericDimensionful.([1,1e-13])
+end
+
 end # module TestGeneric

stdlib/LinearAlgebra/test/givens.jl

@@ -72,32 +72,17 @@ end
 
 # 36430
 # dimensional correctness:
-const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
-isdefined(Main, :Furlongs) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Furlongs.jl"))
-using .Main.Furlongs
 
-@testset "testing dimensions with Furlongs" begin
-    @test_throws MethodError givens(Furlong(1.0), Furlong(2.0), 1, 2)
-end
-
-const TNumber = Union{Float64,ComplexF64}
-struct MockUnitful{T<:TNumber} <: Number
-    data::T
-    MockUnitful(data::T) where T<:TNumber = new{T}(data)
-end
-import Base: *, /, one, oneunit
-*(a::MockUnitful{T}, b::T) where T<:TNumber = MockUnitful(a.data * b)
-*(a::T, b::MockUnitful{T}) where T<:TNumber = MockUnitful(a * b.data)
-*(a::MockUnitful{T}, b::MockUnitful{T}) where T<:TNumber = MockUnitful(a.data * b.data)
-/(a::MockUnitful{T}, b::MockUnitful{T}) where T<:TNumber = a.data / b.data
-one(::Type{<:MockUnitful{T}}) where T = one(T)
-oneunit(::Type{<:MockUnitful{T}}) where T = MockUnitful(one(T))
+@testset "testing dimensions with GenericDimensionfuls" begin
+    @test givens(GenericDimensionful(1.0), GenericDimensionful(2.0), 1, 2) == (givens(1.0,2.0, 1, 2)[1], GenericDimensionful(hypot(1.0,2.0)))
+    @test_throws MethodError givens(GenericDimensionful(1.0), GenericDimensionful{2}(2.0), 1, 2)
 
-@testset "unitful givens rotation unitful $T " for T in (Float64, ComplexF64)
-    g, r = givens(MockUnitful(T(3)), MockUnitful(T(4)), 1, 2)
-    @test g.c ≈ 3/5
-    @test g.s ≈ 4/5
-    @test r.data ≈ 5.0
+    for T in (Float64, ComplexF64)
+        g, r = givens(GenericDimensionful(T(3)), GenericDimensionful(T(4)), 1, 2)
+        @test g.c ≈ 3/5
+        @test g.s ≈ 4/5
+        @test r ≈ GenericDimensionful(5.0)
+    end
 end
 
 end # module TestGivens

stdlib/LinearAlgebra/test/lu.jl

@@ -325,19 +325,16 @@ include("trickyarithmetic.jl")
 end
 
 # dimensional correctness:
-const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
-isdefined(Main, :Furlongs) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Furlongs.jl"))
-using .Main.Furlongs
 
 @testset "lu factorization with dimension type" begin
     n = 4
-    A = Matrix(Furlong(1.0) * I, n, n)
+    A = Matrix(GenericDimensionful(1.0) * I, n, n)
     F = lu(A).factors
     @test Diagonal(F) == Diagonal(A)
-    # upper triangular part has a unit Furlong{1}
-    @test all(x -> typeof(x) == Furlong{1, Float64}, F[i,j] for j=1:n for i=1:j)
-    # lower triangular part is unitless Furlong{0}
-    @test all(x -> typeof(x) == Furlong{0, Float64}, F[i,j] for j=1:n for i=j+1:n)
+    # upper triangular part has a unit GenericDimensionful{1}
+    @test all(x -> typeof(x) == GenericDimensionful{1, Float64}, F[i,j] for j=1:n for i=1:j)
+    # lower triangular part is unitless GenericDimensionful{0}
+    @test all(x -> typeof(x) == GenericDimensionful{0, Float64}, F[i,j] for j=1:n for i=j+1:n)
 end
 
 @testset "Issue #30917. Determinant of integer matrix" begin

stdlib/LinearAlgebra/test/special.jl

@@ -330,9 +330,6 @@ end
 
 
 # for testing types with a dimension
-const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
-isdefined(Main, :Furlongs) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Furlongs.jl"))
-using .Main.Furlongs
 
 @testset "zero and one for structured matrices" begin
     for elty in (Int64, Float64, ComplexF64)
@@ -395,11 +392,12 @@ using .Main.Furlongs
     @test one(S) isa SymTridiagonal
 
     # eltype with dimensions
-    D = Diagonal{Furlong{2, Int64}}([1, 2, 3, 4])
-    Bu = Bidiagonal{Furlong{2, Int64}}([1, 2, 3, 4], [1, 2, 3], 'U')
-    Bl =  Bidiagonal{Furlong{2, Int64}}([1, 2, 3, 4], [1, 2, 3], 'L')
-    T = Tridiagonal{Furlong{2, Int64}}([1, 2, 3], [1, 2, 3, 4], [1, 2, 3])
-    S = SymTridiagonal{Furlong{2, Int64}}([1, 2, 3, 4], [1, 2, 3])
+    GD2 = GenericDimensionful{2}
+    D = Diagonal(GD2.([1, 2, 3, 4]))
+    Bu = Bidiagonal(GD2.([1, 2, 3, 4]), GD2.([1, 2, 3]), 'U')
+    Bl =  Bidiagonal(GD2.([1, 2, 3, 4]), GD2.([1, 2, 3]), 'L')
+    T = Tridiagonal(GD2.([1, 2, 3]), GD2.([1, 2, 3, 4]), GD2.([1, 2, 3]))
+    S = SymTridiagonal(GD2.([1, 2, 3, 4]), GD2.([1, 2, 3]))
     mats = [D, Bu, Bl, T, S]
     for A in mats
         @test iszero(zero(A))

stdlib/LinearAlgebra/test/triangular.jl

@@ -611,13 +611,10 @@ end
 @test UpperTriangular(Matrix(1.0I, 3, 3)) \ view(fill(1., 3), [1,2,3]) == fill(1., 3)
 
 # dimensional correctness:
-const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
-isdefined(Main, :Furlongs) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "Furlongs.jl"))
-using .Main.Furlongs
-LinearAlgebra.sylvester(a::Furlong,b::Furlong,c::Furlong) = -c / (a + b)
+LinearAlgebra.sylvester(a::GenericDimensionful,b::GenericDimensionful,c::GenericDimensionful) = -c / (a + b)
 
-let A = UpperTriangular([Furlong(1) Furlong(4); Furlong(0) Furlong(1)])
-    @test sqrt(A) == Furlong{1//2}.(UpperTriangular([1 2; 0 1]))
+let A = UpperTriangular([GenericDimensionful(1) GenericDimensionful(4); GenericDimensionful(0) GenericDimensionful(1)])
+    @test sqrt(A) == GenericDimensionful{1//2}.(UpperTriangular([1 2; 0 1]))
 end
 
 @testset "similar should preserve underlying storage type" begin

stdlib/Test/docs/src/index.md

@@ -291,3 +291,19 @@ end
 ```@meta
 DocTestSetup = nothing
 ```
+
+## Generic Types
+
+It is often useful to test that packages function correctly with string, container, and number types
+that are similar but not identical to Julia `Base` types like `String`, `Array`, and so on.  The
+`Test` module exports a few such types for testing purposes, so that you don't need to have a test
+dependency on a particular choice of external package implementing an alternative type:
+
+```@docs
+Test.GenericDimensionful
+Test.GenericString
+Test.GenericArray
+Test.GenericDict
+Test.GenericSet
+Test.GenericOrder
+```