Make AbstractDimensions <: Number

[anonymous-shrew]

Background

Hello,

I am the developer of UnitfulTensors.jl, a linear algebra library with support for physical units. At the start of the development, I defined some types for working with scalar quantities, which are basically equivalent to your Dimensions, AbstractDimensions, Quantity, and AbstractQuantity. Now I am considering outsourcing all scalar functionality to DynamicQuantities.jl.

I noticed that things get simpler if AbstractDimensions was a subtype of Number.

Motivation

Physical dimensions are scalar (i. e., non-container) objects supporting arithmetic operations:

• length + length = length, length + time = error
• length - length = length, length - time = error
• length * length = area
• length / length = dimensionless

Other mathematical functions can be defined too:

• abs(length) = length
• exp(dimensionless) = dimensionless, exp(length) = error

With such definitions, most operations on physical quantities can be expressed in a uniform manner:

function f(args::Vararg{Quantity})
    vals = map(ustrip, args)
    dims = map(dimension, args)
    return Quantity(f(vals...), f(dims...))
end

Which is rather neat.

Also, AbstractDimensions <: Number would ensure the correct (scalar-like) broadcasting behavior.

Possible objections

• Addition and subtraction are defined only for identical dimensions.

However, even with more usual Numbers there are cases when some arithmetic operations are undefined. E.g., you cannot divide by zero. You generally cannot divide integers and expect the result to be an integer. For positive integers, even subtraction is generally undefined (if one requires the result to be a positive integer), as well as zero(x). Addition may result in overflow.

• Dimensions cannot be converted to something like Float64 or Complex{Float64}.

But the same can be said about quaternions, transfinite numbers, Grassmann numbers.

And, of course, the same objections would apply to quantities as well.

An unrelated question

Why did you write your own parsing/printing functions instead of using those from Unitful.jl? Defining megameters but not megagrams looks arbitrary. Dimensions are not units, so printing dimension(1u"cm") as "m" is arguably incorrect.

[MilesCranmer]

Interesting! I'm not opposed. Just wondering what this would buy over overloading the individual operators as is done currently? I worry there would be many more functions defined on Number which are not compatible, than functions which are compatible.

[anonymous-shrew]

Doing math with dimensions reduces code repetition, because you can define a lot of operations on quantities in one loop (a real-life example from UnitfulTensors.jl):

for f in (:det, :logdet, :tr, :lyap, :sylvester)
    @eval @inline function $f(args::AbstractUnitfulScalarOrTensor...; kwargs...)
        dims = $f(dimensions.(args)...; kwargs...)
        vals = $f(values.(args)...; kwargs...)
        T = promote_unitful($f, vals, dims, args...)
        return T(vals, dims)
    end
end

and then define only operations on dimensions, without repeating ustrip, dimension, and new_quantity in every function.

And for doing math with dimensions, it is helpful to subtype Number, because then some functions work automatically, including arithmetic with promotion.

I glanced through number.jl in Base and methodswith(Number), and it looks like things that make sense will work, and things that don't make sense will throw a MethodError.

[anonymous-shrew]

Here is an example where I define a simplified version of Dimensions applicable to dimensionless quantities only. Actually, two versions: one is a Number, and the other isn't.

using LinearAlgebra, Test
import Base: +, -, *, /, zero, one, abs, conj, exp, sin, convert, promote_rule

# these structs represent the dimensions of dimensionless quantities
# T represents the type of dimensional exponents
struct DimensionlessAny{T} end
struct DimensionlessNumber{T} <: Number end # despite the name, represents dimensions, not a quantity
const Dimensionless = Union{DimensionlessAny, DimensionlessNumber}

for f in (:+, :-, :*, :/)
    @eval $f(x::T, y::T) where T <: Dimensionless = x
end

for f in (:zero, :abs, :conj, :exp, :sin)
    @eval $f(x::Dimensionless) = x
end

one(::Type{T}) where T <: Dimensionless = T()

# promoting dimensional exponents to a common type
for D in (:DimensionlessAny, :DimensionlessNumber)
    @eval convert(::Type{$D{T1}}, ::$D{T2}) where {T1, T2} = $D{T1}()
    @eval promote_rule(::Type{$D{T1}}, ::Type{$D{T2}}) where {T1, T2} = $D{promote_type(T1, T2)}
end

# silence errors
Test.record(ts::Test.DefaultTestSet, t::Union{Test.Error, Test.Fail}) = (push!(ts.results, t); t)

@testset "all" begin
    for (desc, D) in (("Any"   , DimensionlessAny   ),
                      ("Number", DimensionlessNumber))
        nodims = D{Float64}(); nodimsInt = D{Int}(); T = D{Float64}
        @testset "$desc" begin
            # misc
            @test T(nodims) == nodims
            @test copy(nodims) == nodims
            # math
            @test +(nodims) == nodims
            @test *(nodims) == nodims
            @test inv(nodims) == nodims
            @test nodims^2 == nodims
            @test sign(nodims) == nodims
            @test abs2(nodims) == nodims
            @test ℯ^nodims == nodims
            @test csc(nodims) == nodims
            @test nodims + missing === missing
            @test nodims * [nodims] == [nodims]
            # promotion
            @test nodims == nodimsInt
            @test nodims + nodimsInt == nodims
            @test nodims - nodimsInt == nodims
            @test nodims * nodimsInt == nodims
            @test nodims / nodimsInt == nodims
            # array-like behavior
            @test eltype(nodims) == T
            @test ndims(nodims) == 0
            @test size(nodims) == ()
            @test axes(nodims) == ()
            @test length(nodims) == 1
            @test firstindex(nodims) == 1
            @test first(nodims) == nodims
            @test nodims[] == nodims
            @test get(nodims, (1, 1), nothing) == nodims
            # collection-like behavior
            @test Base.IteratorSize(T) == Base.HasShape{0}()
            @test iterate(nodims) == (nodims, nothing)
            @test iterate(nodims, nothing) == nothing
            @test !isempty(nodims)
            @test nodims ∈ nodims
            @test map(+, nodims, nodims) == nodims + nodims
            # matrix-like behavior
            @test transpose(nodims) == nodims
            @test adjoint(nodims) == nodims
            @test det(nodims) == nodims
            @test tr(nodims) == nodims
            @test issymmetric(nodims) == true
            @test LinearAlgebra.symmetric(nodims) == nodims
            # broadcasting
            @test nodims .* [nodims] == [nodims]
        end
    end
end

Everything works for DimensionlessNumber, but nothing for DimensionlessAny.

Calculating the determinant of a number may not be very useful, but it's in stdlib for some reason, so why not make it work with AbstractDimensions.

[anonymous-shrew]

And here is an example where functions that are not applicable to dimensions just throw MethodErrors:

nodims = DimensionlessNumber{Float64}()
isodd(nodims)
# MethodError: no method matching isreal(::DimensionlessNumber{Float64})
isposdef(nodims)
# MethodError: no method matching imag(::DimensionlessNumber{Float64})
AbstractChar(nodims)
# MethodError: no method matching UInt32(::DimensionlessNumber{Float64})