Create in-house units module

[MilesCranmer]

This adds an in-house units module, so now you can parse a wide variety of base and derived SI units without relying on Unitful.jl. It includes a few non SI units: min, h, day, year - but I purposefully didn't go too far beyond that, to avoid potential namespace ambiguities. I did not include any Gaussian units (like erg - which is used in astronomy), because Gaussian units require different formulas and I feel like mixing the unit systems might lead to confusion.

This PR also makes it so that all physical quantities are displayed in SI base units "m¹ᐟ² kg¹ mol³", rather than the unevaluated "𝐋 ¹ᐟ² 𝐌 ¹ 𝐍 ³" representation, so representation is more obvious.

Here are all the available units and their prefixes:

# SI base units
const m = Quantity(1.0, length=1)
const g = Quantity(1e-3, mass=1)
const s = Quantity(1.0, time=1)
const A = Quantity(1.0, current=1)
const K = Quantity(1.0, temperature=1)
const cd = Quantity(1.0, luminosity=1)
const mol = Quantity(1.0, amount=1)

@add_prefixes m (f, p, n, μ, u, c, m, k, M, G)
@add_prefixes g (μ, u, m, k)
@add_prefixes s (f, p, n, μ, u, m)
@add_prefixes A (n, μ, u, m, k)
@add_prefixes K (m,)
@add_prefixes cd (m,)
@add_prefixes mol (m,)

# SI derived units
const Hz = inv(s)
const N = kg * m / s^2
const Pa = N / m^2
const J = N * m
const W = J / s
const C = A * s
const V = W / A
const F = C / V
const Ω = V / A
const T = N / (A * m)

@add_prefixes Hz (k, M, G)
@add_prefixes N ()
@add_prefixes P ()
@add_prefixes J ()
@add_prefixes W (k, M, G)
@add_prefixes C ()
@add_prefixes V (m, k, M, G)
@add_prefixes F ()
@add_prefixes Ω (m,)
@add_prefixes T ()

# Assorted units
const min = 60 * s
const h = 60 * min
const day = 24 * h
const yr = 365.25 * day

I tried to add all of the common units. For physical constants or other weird unit systems, I think it is best we make the user define them themselves, to ensure they know what they are doing, and avoid potential namespace collisions.

You can access these units with uparse or @u_str, as in Unitful. For example:

julia> x = 0.5u"km/s"
500.0 m s⁻¹

julia> x = 100u"ug/day"
1.1574074074074076e-12 kg s⁻¹

julia> x = 100u"fm^2"
1.0000000000000001e-28 m²

julia> x = 0.5u"kW*h"
1.8e6 m² kg s⁻²

You can still convert from Unitful.jl if you need to, but I think this is a more robust option.

What do you think? @oscardssmith @j-fu @odow

---

TODO:

• Throw an error if the user is converting from Unitful.jl with a non-SI upreferred.
• Update README
• Add docs on available units

[github-actions]

Benchmark Results

	main	37de814...	t[main]/t[37de814...]
creation/Quantity(x)	3.4 ± 0 ns	3.4 ± 0 ns	1
creation/Quantity(x, length=y)	3.4 ± 0 ns	3.4 ± 0 ns	1
time_to_load	0.0849 ± 0.0021 s	0.101 ± 0.0015 s	0.841
with_numbers/*real	3.1 ± 0.1 ns	3.4 ± 0 ns	0.912
with_numbers/^int	0.0331 ± 0.0025 μs	23.8 ± 2.4 ns	1.39
with_numbers/^int * real	0.0344 ± 0.0025 μs	23.8 ± 2.4 ns	1.45
with_quantity/+y	4.3 ± 0 ns	7.7 ± 0.1 ns	0.558
with_quantity//y	1.7 ± 0 ns	4.3 ± 0.1 ns	0.395
with_self/dimension	1.7 ± 0 ns	1.7 ± 0 ns	1
with_self/inv	1.7 ± 0 ns	3.7 ± 0.1 ns	0.459
with_self/ustrip	1.7 ± 0 ns	1.7 ± 0 ns	1

Benchmark Plots

A plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR.
Go to "Actions"->"Benchmark a pull request"->[the most recent run]->"Artifacts" (at the bottom).

[MilesCranmer]

I'm debating whether to add c, the speed of light, as a derived unit. I'm hesitant because c feels more like a physical constant, and I wouldn't want to mislead people that physical constants are present: if c is included, people might think h (planck constant) is also included. But in actuality, h = hours.

I think I won't include it, but just wanted to hear if others have thoughts.

The other non-SI units I added are L (liters), bar, and eV - all of which are fairly common in scientific calculations. eV is technically a "constant" but due to its usage it feels more like a unit. (eV, year, and day are probably the only ones which can get away with being "units" despite technically being physical constants.)

[MilesCranmer]

Quick update: after some good points made by @mcabbott, this PR now sets units to the type Quantity{Rational{Int64},...}. This gives you better promotion properties, like to Float16, whereas before it would automatically convert everything to Float64:

julia> typeof(u"1")
Quantity{Int64, DynamicQuantities.FixedRational{Int32, 25200}}

julia> typeof(u"km")
Quantity{Rational{Int64}, DynamicQuantities.FixedRational{Int32, 25200}}

julia> typeof(Float16(1.0)u"km")
Quantity{Float16, DynamicQuantities.FixedRational{Int32, 25200}}

julia> typeof(1f0u"km")
Quantity{Float32, DynamicQuantities.FixedRational{Int32, 25200}}

julia> typeof(1.0u"km")
Quantity{Float64, DynamicQuantities.FixedRational{Int32, 25200}}

---

Edit: the one downside is that we can't store meV as a unit. So to make things simple I'm just going to take out eV altogether as its more of a physical constant.

[MilesCranmer]

Hm, the one issue with this is you can get overflow errors for certain units:

julia> u"fm^2"
ERROR: LoadError: OverflowError: 1000000000000000 * 1000000000000000 overflowed for type Int64

$\text{fm}^2$ is actually used for scattering cross sections in particle physics: https://en.wikipedia.org/wiki/Barn_(unit)

@mcabbott

[MilesCranmer]

Reverted units back to Float64.

[MilesCranmer]

I think we can merge this in a couple of days; let me know if there is any other feedback.
Cheers,
Miles

[MilesCranmer]

@mcabbott Maybe something like this could be an option for the type of each unit? It's a Float64 that tries to demote itself whenever possible:

struct WeakFloat <: AbstractFloat
    x::Float64
end

Base.promote(x::WeakFloat, y::Float64) = (x.x, y)
Base.promote(x::WeakFloat, y::AbstractFloat) = (convert(typeof(y), x.x), y)
Base.promote(x::WeakFloat, y) = promote(x.x, y)
Base.promote(x, y::WeakFloat) = reverse(promote(y, x))

function Base.promote(x::WeakFloat, y::Rational{R}) where {R}
    # If under eps * 1e5, we promote to a float:
    abs(x.x) < eps(Float64) * 1e5 && return promote(x.x, y)
    return (rationalize(R, x.x), y)
end

function Base.promote(x::WeakFloat, y::Integer)
    abs(x.x) < eps(Float64) * 1e5 && return promote(x.x, y)
    let r=rationalize(typeof(y), x.x)
        return isinteger(r) ? (convert(typeof(y), r), y) : promote(r, y)
    end
end

Which gives you:

julia> 2 * WeakFloat(100.0)
200

julia> (2//1) * WeakFloat(0.5)
1//1

julia> WeakFloat(1e-5) * 10000
1//10

julia> WeakFloat(1e-20) * 20
2.0e-19

The downside is that this is type unstable: the promoted type varies depending on how close the float is to 0.

[MilesCranmer]

Thanks for the reviews, everyone!