SymbolicML · MilesCranmer · Jun 29, 2023 · Jun 29, 2023 · Jul 9, 2023 · Jul 9, 2023
diff --git a/src/DynamicQuantities.jl b/src/DynamicQuantities.jl
@@ -8,6 +8,8 @@ export ulength, umass, utime, ucurrent, utemperature, uluminosity, uamount
 export uparse, @u_str, sym_uparse, @us_str, uexpand, uconvert
 
 include("fixed_rational.jl")
+include("lazy_float.jl")
+
 include("types.jl")
 include("utils.jl")
 include("math.jl")
@@ -22,6 +24,7 @@ export expand_units
 
 import PackageExtensionCompat: @require_extensions
 import .Units
+import .Units: DEFAULT_UNIT_TYPE
 import .Constants
 import .UnitsParse: uparse, @u_str
 

diff --git a/src/constants.jl b/src/constants.jl
@@ -1,12 +1,11 @@
 module Constants
 
-import ..DEFAULT_QUANTITY_TYPE
 import ..Quantity
 import ..Units as U
-import ..Units: _add_prefixes
+import ..Units: _add_prefixes, DEFAULT_UNIT_TYPE
 
 const _CONSTANT_SYMBOLS = Symbol[]
-const _CONSTANT_VALUES = DEFAULT_QUANTITY_TYPE[]
+const _CONSTANT_VALUES = DEFAULT_UNIT_TYPE[]
 
 macro register_constant(name, value)
     return esc(_register_constant(name, value))
@@ -20,7 +19,7 @@ end
 function _register_constant(name::Symbol, value)
     s = string(name)
     return quote
-        const $name = $value
+        const $name = convert(DEFAULT_UNIT_TYPE, $value)
         push!(_CONSTANT_SYMBOLS, Symbol($s))
         push!(_CONSTANT_VALUES, $name)
     end
@@ -87,7 +86,7 @@ end
 )
 
 # Measured
-@register_constant alpha DEFAULT_QUANTITY_TYPE(7.2973525693e-3)
+@register_constant alpha DEFAULT_UNIT_TYPE(7.2973525693e-3)
 @register_constant u 1.66053906660e-27 * U.kg
 @register_constant G 6.67430e-11 * U.m^3 / (U.kg * U.s^2)
 @register_constant mu_0 4π * alpha * hbar / (e^2 * c)

diff --git a/src/lazy_float.jl b/src/lazy_float.jl
@@ -0,0 +1,43 @@
+"""
+    AutoFloat <: AbstractFloat
+
+A wrapper around a `Float64` which automatically demotes
+itself to `Float32` or `Float16` if interacting with
+such a type. In other cases it will respect promotion
+rules associated with `Float64`.
+"""
+struct AutoFloat <: AbstractFloat
+    value::Float64
+
+    AutoFloat(x::AbstractFloat) = new(convert(Float64, x))
+end
+
+Base.float(x::AutoFloat) = x.value
+
+Base.convert(::Type{AutoFloat}, x::AutoFloat) = x
+Base.convert(::Type{AutoFloat}, x::FixedRational) = AutoFloat(convert(Float64, x))
+Base.convert(::Type{AutoFloat}, x::Number) = AutoFloat(x)
+Base.convert(::Type{T}, x::AutoFloat) where {T<:Number} = convert(T, float(x))
+Base.promote_rule(::Type{AutoFloat}, ::Type{T}) where {T<:AbstractFloat} = T
+Base.promote_rule(::Type{AutoFloat}, ::Type{T}) where {T} = promote_type(Float64, T)
+
+Base.show(io::IO, x::AutoFloat) = print(io, float(x))
+
+Base.:(==)(a::AutoFloat, b::AutoFloat) = float(a) == float(b)
+Base.:+(a::AutoFloat, b::AutoFloat) = AutoFloat(float(a) + float(b))
+Base.:-(a::AutoFloat) = AutoFloat(-float(a))
+Base.:-(a::AutoFloat, b::AutoFloat) = AutoFloat(float(a) + float(-b))
+Base.:*(a::AutoFloat, b::AutoFloat) = AutoFloat(float(a) * float(b))
+Base.inv(a::AutoFloat) = AutoFloat(inv(float(a)))
+Base.abs(a::AutoFloat) = AutoFloat(abs(float(a)))
+Base.:/(a::AutoFloat, b::AutoFloat) = a * inv(b)
+Base.:^(a::AutoFloat, b::Int) = AutoFloat(float(a) ^ b)
+Base.:^(a::AutoFloat, b::AutoFloat) = AutoFloat(float(a) ^ float(b))
+Base.sqrt(a::AutoFloat) = AutoFloat(sqrt(float(a)))
+Base.cbrt(a::AutoFloat) = AutoFloat(cbrt(float(a)))
+Base.eps(::Type{AutoFloat}) = eps(Float64)
+
+# Ambiguities:
+for T in (:(Rational{<:Any}), :(Base.TwicePrecision), :AbstractChar, :Complex, :Number)
+    @eval AutoFloat(x::$T) = AutoFloat(float(x))
+end
diff --git a/src/symbolic_dimensions.jl b/src/symbolic_dimensions.jl
@@ -1,5 +1,6 @@
-import .Units: UNIT_SYMBOLS, UNIT_MAPPING, UNIT_VALUES
+import .Units: UNIT_SYMBOLS, UNIT_MAPPING, UNIT_VALUES, DEFAULT_UNIT_BASE_TYPE
 import .Constants: CONSTANT_SYMBOLS, CONSTANT_MAPPING, CONSTANT_VALUES
+import ..DEFAULT_DIM_BASE_TYPE
 
 const SYMBOL_CONFLICTS = intersect(UNIT_SYMBOLS, CONSTANT_SYMBOLS)
 
@@ -93,6 +94,8 @@ function Base.convert(::Type{Quantity{T,D}}, q::Quantity{<:Any,<:SymbolicDimensi
     return result
 end
 
+const DEFAULT_SYMBOLIC_UNIT_TYPE = Quantity{DEFAULT_UNIT_BASE_TYPE,SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}}
+
 """
     uexpand(q::Quantity{<:Any,<:SymbolicDimensions})
 
@@ -269,20 +272,22 @@ module SymbolicUnitsParse
     import ..CONSTANT_SYMBOLS
     import ..SYMBOL_CONFLICTS
     import ..SymbolicDimensions
+    import ..DEFAULT_SYMBOLIC_UNIT_TYPE
+    import ..DEFAULT_UNIT_BASE_TYPE
+    import ..DEFAULT_DIM_BASE_TYPE
 
     import ...Quantity
-    import ...DEFAULT_VALUE_TYPE
-    import ...DEFAULT_DIM_BASE_TYPE
 
     # Lazily create unit symbols (since there are so many)
     module Constants
         import ..CONSTANT_SYMBOLS
         import ..SYMBOL_CONFLICTS
         import ..SymbolicDimensions
+        import ..DEFAULT_SYMBOLIC_UNIT_TYPE
+        import ..DEFAULT_UNIT_BASE_TYPE
+        import ..DEFAULT_DIM_BASE_TYPE
 
         import ..Quantity
-        import ..DEFAULT_VALUE_TYPE
-        import ..DEFAULT_DIM_BASE_TYPE
 
         import ...Constants as EagerConstants
 
@@ -292,11 +297,11 @@ module SymbolicUnitsParse
             CONSTANT_SYMBOLS_EXIST[] || lock(CONSTANT_SYMBOLS_LOCK) do
                 CONSTANT_SYMBOLS_EXIST[] && return nothing
                 for unit in setdiff(CONSTANT_SYMBOLS, SYMBOL_CONFLICTS)
-                    @eval const $unit = Quantity(DEFAULT_VALUE_TYPE(1.0), SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}; $(unit)=1)
+                    @eval const $unit = Quantity(DEFAULT_UNIT_BASE_TYPE(1.0), SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}; $(unit)=1)
                 end
                 # Evaluate conflicting symbols to non-symbolic form:
                 for unit in SYMBOL_CONFLICTS
-                    @eval const $unit = convert(Quantity{DEFAULT_VALUE_TYPE,SymbolicDimensions}, EagerConstants.$unit)
+                    @eval const $unit = convert(DEFAULT_SYMBOLIC_UNIT_TYPE, EagerConstants.$unit)
                 end
                 CONSTANT_SYMBOLS_EXIST[] = true
             end
@@ -311,7 +316,7 @@ module SymbolicUnitsParse
         UNIT_SYMBOLS_EXIST[] || lock(UNIT_SYMBOLS_LOCK) do
             UNIT_SYMBOLS_EXIST[] && return nothing
             for unit in UNIT_SYMBOLS
-                @eval const $unit = Quantity(DEFAULT_VALUE_TYPE(1.0), SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}; $(unit)=1)
+                @eval const $unit = Quantity(DEFAULT_UNIT_BASE_TYPE(1.0), SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}; $(unit)=1)
             end
             UNIT_SYMBOLS_EXIST[] = true
         end
@@ -338,11 +343,11 @@ module SymbolicUnitsParse
         _generate_unit_symbols()
         Constants._generate_unit_symbols()
         raw_result = eval(Meta.parse(raw_string))
-        return copy(as_quantity(raw_result))::Quantity{DEFAULT_VALUE_TYPE,SymbolicDimensions{DEFAULT_DIM_BASE_TYPE}}
+        return copy(as_quantity(raw_result))::DEFAULT_SYMBOLIC_UNIT_TYPE
     end
 
-    as_quantity(q::Quantity) = q
-    as_quantity(x::Number) = Quantity(convert(DEFAULT_VALUE_TYPE, x), SymbolicDimensions{DEFAULT_DIM_BASE_TYPE})
+    as_quantity(q::Quantity) = convert(DEFAULT_SYMBOLIC_UNIT_TYPE, q)
+    as_quantity(x::Number) = DEFAULT_SYMBOLIC_UNIT_TYPE(x)
     as_quantity(x) = error("Unexpected type evaluated: $(typeof(x))")
 end
 

diff --git a/src/units.jl b/src/units.jl
@@ -2,13 +2,15 @@ module Units
 
 import ..DEFAULT_DIM_TYPE
 import ..DEFAULT_VALUE_TYPE
-import ..DEFAULT_QUANTITY_TYPE
 import ..Quantity
+import ..Quantity
+import ..AutoFloat
 
-@assert DEFAULT_VALUE_TYPE == Float64 "`units.jl` must be updated to support a different default value type."
+const DEFAULT_UNIT_BASE_TYPE = AutoFloat
+const DEFAULT_UNIT_TYPE = Quantity{DEFAULT_UNIT_BASE_TYPE,DEFAULT_DIM_TYPE}
 
 const _UNIT_SYMBOLS = Symbol[]
-const _UNIT_VALUES = DEFAULT_QUANTITY_TYPE[]
+const _UNIT_VALUES = DEFAULT_UNIT_TYPE[]
 
 macro register_unit(name, value)
     return esc(_register_unit(name, value))
@@ -22,7 +24,7 @@ end
 function _register_unit(name::Symbol, value)
     s = string(name)
     return quote
-        const $name = $value
+        const $name = convert(DEFAULT_UNIT_TYPE, $value)
         push!(_UNIT_SYMBOLS, Symbol($s))
         push!(_UNIT_VALUES, $name)
     end
@@ -37,19 +39,19 @@ function _add_prefixes(base_unit::Symbol, prefixes, register_function)
     for (prefix, value) in zip(keys(all_prefixes), values(all_prefixes))
         prefix in prefixes || continue
         new_unit = Symbol(prefix, base_unit)
-        push!(expr.args, register_function(new_unit, :($value * $base_unit)))
+        push!(expr.args, register_function(new_unit, :(convert(DEFAULT_UNIT_TYPE, $value * $base_unit))))
     end
     return expr
 end
 
 # SI base units
-@register_unit m DEFAULT_QUANTITY_TYPE(1.0, length=1)
-@register_unit g DEFAULT_QUANTITY_TYPE(1e-3, mass=1)
-@register_unit s DEFAULT_QUANTITY_TYPE(1.0, time=1)
-@register_unit A DEFAULT_QUANTITY_TYPE(1.0, current=1)
-@register_unit K DEFAULT_QUANTITY_TYPE(1.0, temperature=1)
-@register_unit cd DEFAULT_QUANTITY_TYPE(1.0, luminosity=1)
-@register_unit mol DEFAULT_QUANTITY_TYPE(1.0, amount=1)
+@register_unit m Quantity(1.0, length=1)
+@register_unit g Quantity(1e-3, mass=1)
+@register_unit s Quantity(1.0, time=1)
+@register_unit A Quantity(1.0, current=1)
+@register_unit K Quantity(1.0, temperature=1)
+@register_unit cd Quantity(1.0, luminosity=1)
+@register_unit mol Quantity(1.0, amount=1)
 
 @add_prefixes m (f, p, n, μ, u, c, d, m, k, M, G)
 @add_prefixes g (μ, u, m, k)

diff --git a/src/uparse.jl b/src/uparse.jl
@@ -3,7 +3,7 @@ module UnitsParse
 import ..Quantity
 import ..DEFAULT_DIM_TYPE
 import ..DEFAULT_VALUE_TYPE
-import ..Units: UNIT_SYMBOLS
+import ..Units: UNIT_SYMBOLS, DEFAULT_UNIT_TYPE
 import ..Constants
 
 function _generate_units_import()
@@ -30,11 +30,11 @@ the quantity corresponding to the speed of light multiplied by Hertz,
 squared.
 """
 function uparse(s::AbstractString)
-    return as_quantity(eval(Meta.parse(s)))::Quantity{DEFAULT_VALUE_TYPE,DEFAULT_DIM_TYPE}
+    return as_quantity(eval(Meta.parse(s)))::DEFAULT_UNIT_TYPE
 end
 
-as_quantity(q::Quantity) = q
-as_quantity(x::Number) = Quantity(convert(DEFAULT_VALUE_TYPE, x), DEFAULT_DIM_TYPE)
+as_quantity(q::Quantity) = convert(DEFAULT_UNIT_TYPE, q)
+as_quantity(x::Number) = DEFAULT_UNIT_TYPE(x)
 as_quantity(x) = error("Unexpected type evaluated: $(typeof(x))")
 
 """

diff --git a/test/unittests.jl b/test/unittests.jl
@@ -1,13 +1,20 @@
 using DynamicQuantities
 using DynamicQuantities: FixedRational
-using DynamicQuantities: DEFAULT_DIM_BASE_TYPE, DEFAULT_DIM_TYPE, DEFAULT_VALUE_TYPE
+using DynamicQuantities: DEFAULT_DIM_BASE_TYPE, DEFAULT_DIM_TYPE, DEFAULT_VALUE_TYPE, DEFAULT_UNIT_TYPE
+using DynamicQuantities: AutoFloat
 using DynamicQuantities: array_type, value_type, dim_type, quantity_type
 using Ratios: SimpleRatio
 using SaferIntegers: SafeInt16
 using StaticArrays: SArray, MArray
 using LinearAlgebra: norm
 using Test
 
+function show_string(i)
+    io = IOBuffer()
+    show(io, i)
+    return String(take!(io))
+end
+
 @testset "Basic utilities" begin
 
     for T in [DEFAULT_VALUE_TYPE, Float16, Float32, Float64], R in [DEFAULT_DIM_BASE_TYPE, Rational{Int16}, Rational{Int32}, SimpleRatio{Int}, SimpleRatio{SafeInt16}]
@@ -379,17 +386,52 @@ end
     @test ustrip(z) ≈ 60 * 60 * 24 * 365.25
 
     # Test type stability of extreme range of units
-    @test typeof(u"1") == Quantity{Float64,DEFAULT_DIM_TYPE}
-    @test typeof(u"1f0") == Quantity{Float64,DEFAULT_DIM_TYPE}
-    @test typeof(u"s"^2) == Quantity{Float64,DEFAULT_DIM_TYPE}
-    @test typeof(u"Ω") == Quantity{Float64,DEFAULT_DIM_TYPE}
-    @test typeof(u"Gyr") == Quantity{Float64,DEFAULT_DIM_TYPE}
-    @test typeof(u"fm") == Quantity{Float64,DEFAULT_DIM_TYPE}
-    @test typeof(u"fm"^2) == Quantity{Float64,DEFAULT_DIM_TYPE}
+    @test typeof(u"1") == DEFAULT_UNIT_TYPE
+    @test typeof(u"1f0") == DEFAULT_UNIT_TYPE
+    @test typeof(u"s"^2) == DEFAULT_UNIT_TYPE
+    @test typeof(u"Ω") == DEFAULT_UNIT_TYPE
+    @test typeof(u"Gyr") == DEFAULT_UNIT_TYPE
+    @test typeof(u"fm") == DEFAULT_UNIT_TYPE
+    @test typeof(u"fm"^2) == DEFAULT_UNIT_TYPE
+
+    # Test type demotion
+    @test typeof(1u"m") == Quantity{Float64,DEFAULT_DIM_TYPE}
+    @test typeof(1f0u"m") == Quantity{Float32,DEFAULT_DIM_TYPE}
+    @test typeof(1.0u"m") == Quantity{Float64,DEFAULT_DIM_TYPE}
+    @test typeof(Float16(1.0)u"m") == Quantity{Float16,DEFAULT_DIM_TYPE}
+
+    @test typeof(1u"m^2/s") == Quantity{Float64,DEFAULT_DIM_TYPE}
+    @test typeof(1f0u"m^2/s") == Quantity{Float32,DEFAULT_DIM_TYPE}
+    @test typeof(1.0u"m^2/s") == Quantity{Float64,DEFAULT_DIM_TYPE}
 
     @test_throws LoadError eval(:(u":x"))
 end
 
+@testset "AutoFloat" begin
+    @test promote_type(AutoFloat, Float16) == Float16
+    @test promote_type(AutoFloat, Float32) == Float32
+    @test promote_type(AutoFloat, Float64) == Float64
+    @test promote_type(AutoFloat, BigFloat) == BigFloat
+    @test promote_type(AutoFloat, Int64) == Float64
+    @test promote_type(AutoFloat, ComplexF16) == promote_type(Float64, ComplexF16)
+
+    x = AutoFloat(1.5)
+    @test show_string(x) == "1.5"
+
+    @test -x == AutoFloat(-1.5)
+    @test abs(-x) == x
+    @test sqrt(x) == AutoFloat(sqrt(1.5))
+    @test cbrt(x) == AutoFloat(cbrt(1.5))
+    @test inv(x) == AutoFloat(inv(1.5))
+
+    y = AutoFloat(2.1)
+    @test x + y == AutoFloat(1.5 + 2.1)
+    @test x - y == AutoFloat(1.5 - 2.1)
+
+    # Should promote to array:
+    @test typeof([u"km/s", 1.5u"km/s"]) <: Vector{<:Quantity{Float64}}
+end
+
 @testset "Constants" begin
     @test Constants.h * Constants.c / (1000.0u"nm") ≈ 1.9864458571489284e-19u"J"
 
@@ -421,11 +463,6 @@ end
     @test convert(Rational, FixedRational{UInt8,6}(2)) === Rational{UInt8}(2)
 
     # Showing rationals
-    function show_string(i)
-        io = IOBuffer()
-        show(io, i)
-        return String(take!(io))
-    end
     @test show_string(FixedRational{Int,10}(2)) == "2"
     @test show_string(FixedRational{Int,10}(11//10)) == "11//10"
 
@@ -642,7 +679,7 @@ end
     @test promote(x, y) == (x, y)
     @test_throws ErrorException promote(x, convert(FixedRational{Int32,100}, 10))
     @test round(Missing, x) === missing
-    @test promote_type(typeof(u"km/s"), typeof(convert(Quantity{Float32}, u"km/s"))) <: Quantity{Float64}
+    @test promote_type(typeof(1.0u"km/s"), typeof(convert(Quantity{Float32}, u"km/s"))) <: Quantity{Float64}
 
     x = 1.0u"m"
     y = missing