Unitful.jl without parametrizing types on units?

General Usage
,
1

I am trying to add support for dimensional constraints in SymbolicRegression.jl: https://github.com/MilesCranmer/SymbolicRegression.jl/pull/220. However, I realize now this is going to be impossible with Unitful.jl due to the fact that each combination of units will have a different type.

Does anybody know if there is a version of Unitful.jl (or alternative units package) that does not parametrize quantity types by units? In other words, I want the physical quantity

1.0u"m/s"

(one meter per second) to be of type

Quantity{Float64}

rather than type:

Quantity{Float64, 𝐋 𝐓⁻¹, Unitful.FreeUnits{(m, s⁻¹), 𝐋 𝐓⁻¹, nothing}}

I want quantities to store their units and dimensions in a statically-typed field, rather than type parameters.

The reason for this is because SymbolicRegression.jl searches over arbitrary expressions, like so:

This means that there are arbitrary combinations of units that we can encounter! So the number of compiled methods will explode and type inference will bottleneck the evaluation. Only with statically-typed quantities will this be able to work.

3 Likes
2

There isn’t one, but I’ve been saying for a few years that someone needs to make such a library. So it’s a nice open potential project.

1 Like
3

It would be nice if a new project keeps the same interface as Unitful as much as possible, so users can easily switch between them.

4 Likes
4

Okay I messed around a bit with getting something barebones to work, just so I can do dimensional analysis… Here’s my attempt. It uses a Dict{Symbol,Rational{Int}} to store the dimensions. It doesn’t know about units, but maybe one could use Unitful as the interface, convert to SI units, strip the units, and then do manipulations with something like this?

const D_TYPE = Dict{Symbol,Rational{Int}}
const COMPAT_D_TYPE = Union{D_TYPE,Vector{Pair{Symbol,Rational{Int}}},Vector{Pair{Symbol,Int}}}
struct Dimensions
    data::D_TYPE

    Dimensions() = new(D_TYPE())
    Dimensions(data::D_TYPE) = new(data)
    Dimensions(data::COMPAT_D_TYPE) = new(D_TYPE(data))
end
struct Quantity{T}
    val::T
    dimensions::Dimensions
    valid::Bool

    Quantity(x) = new{typeof(x)}(x, Dimensions(), true)
    Quantity(x, dimensions::Dimensions) = new{typeof(x)}(x, dimensions, true)
    Quantity(x, data::COMPAT_D_TYPE) = new{typeof(x)}(x, Dimensions(data), true)
    Quantity(x, valid::Bool) = new{typeof(x)}(x, Dimensions(), valid)
    Quantity(x, dimensions::Dimensions, valid::Bool) = new{typeof(x)}(x, dimensions, valid)
    Quantity(x, data::COMPAT_D_TYPE, valid::Bool) = new{typeof(x)}(x, Dimensions(data), valid)
end

Base.show(io::IO, d::Dimensions) = foreach(k -> d[k] != 0 ? (print(io, k); pretty_print_exponent(io, d[k]); print(io, " ")) : print(io, ""), keys(d))
Base.show(io::IO, q::Quantity) = q.valid ? print(io, q.val, " ", q.dimensions) : print(io, "INVALID")
tryround(x::Rational{Int}) = isinteger(x) ? round(Int, x) : x
pretty_print_exponent(io::IO, x::Rational{Int}) = (x >= 0 && isinteger(x)) ? print(io, "^", round(Int, x)) : print(io, "^(", tryround(x), ")")
Base.convert(::Type{Dimensions}, d::D_TYPE) = Dimensions(d)
Base.isfinite(q::Quantity) = isfinite(q.val)
Base.keys(d::Dimensions) = keys(d.data)
Base.values(d::Dimensions) = values(d.data)
Base.iszero(d::Dimensions) = all(iszero, values(d))
Base.getindex(d::Dimensions, k::Symbol) = get(d.data, k, zero(Rational{Int}))
Base.:(==)(l::Dimensions, r::Dimensions) = all(k -> (l[k] == r[k]), union(keys(l), keys(r)))
Base.convert(::Type{T}, q::Quantity) where {T<:Real} = q.valid ? (iszero(q.dimensions) ? convert(T, q.val) : throw(error("Quantity $(q) has dimensions!"))) : throw(error("Quantity $(q) is invalid!"))
Base.float(q::Quantity{T}) where {T<:AbstractFloat} = convert(T, q)

Base.:*(l::Dimensions, r::Dimensions) = Dimensions(D_TYPE([(k, l[k] + r[k]) for k in union(keys(l.data), keys(r.data))]))
Base.:/(l::Dimensions, r::Dimensions) = Dimensions(D_TYPE([(k, l[k] - r[k]) for k in union(keys(l.data), keys(r.data))]))
Base.inv(d::Dimensions) = Dimensions(D_TYPE([(k, -d[k]) for k in keys(d.data)]))

Base.:*(l::Quantity, r::Quantity) = Quantity(l.val * r.val, l.dimensions * r.dimensions, l.valid && r.valid)
Base.:/(l::Quantity, r::Quantity) = Quantity(l.val / r.val, l.dimensions / r.dimensions, l.valid && r.valid)
Base.:+(l::Quantity, r::Quantity) = Quantity(l.val + r.val, l.dimensions, l.dimensions == r.dimensions)
Base.:-(l::Quantity, r::Quantity) = Quantity(l.val - r.val, l.dimensions, l.dimensions == r.dimensions)
Base.:^(l::Quantity, r::Quantity) = let rr=convert(Rational{Int}, r.val); Quantity(l.val ^ rr, l.dimensions ^ rr, l.valid && r.valid && iszero(r.dimensions)); end

Base.:*(l::Quantity, r::Number) = Quantity(l.val * r, l.dimensions, l.valid)
Base.:*(l::Number, r::Quantity) = Quantity(l * r.val, r.dimensions, r.valid)
Base.:/(l::Quantity, r::Number) = Quantity(l.val / r, l.dimensions, l.valid)
Base.:/(l::Number, r::Quantity) = l * inv(r)
Base.:^(l::Dimensions, r::Rational{Int}) = Dimensions(D_TYPE([(k, l.data[k] * r) for k in keys(l.data)]))
Base.:^(l::Quantity, r::Number) = let rr=convert(Rational{Int}, r); Quantity(l.val ^ rr, l.dimensions ^ rr, l.valid); end
Base.inv(q::Quantity) = Quantity(inv(q.val), inv(q.dimensions), q.valid)
Base.sqrt(q::Quantity) = Quantity(sqrt(q.val), q.dimensions ^ (1//2), q.valid)

This also avoids having to throw DimensionError by just having the .valid field. Then you can do stuff like:

julia> x = Quantity(0.2, [:M => 1, :L => 1//2])
0.2 M^1 L^(1//2) 

julia> y = Quantity(10.2, [:M => 2, :T => -2])
10.2 T^(-2) M^2 

julia> x * y
2.04 T^(-2) M^3 L^(1//2)

julia> x / y
0.019607843137254905 T^2 M^(-1) L^(1//2)

julia> x ^ 3
0.008000000000000002 M^3 L^(3//2) 

julia> x ^ -1
5.0 M^(-1) L^(-1//2)

julia> sqrt(x)
0.4472135954999579 M^(1//2) L^(1//4)

julia> x ^ 1.5
0.0894427190999916 M^(3//2) L^(3//4)

and check for dimensionality issues with

julia> x + 2 * x
0.6000000000000001 M^1 L^(1//2)

julia> x + y
INVALID

we can see the second one has .valid==false (but doesn’t throw an expensive error).

Finally, we can have type-stable vectors with different units:

julia> v = [Quantity(randn(), [:M=>rand(0:10)//rand(1:10), :L=>rand(0:10)//rand(1:10)]) for _=1:5]
5-element Vector{Quantity{Float64}}:
 0.8249581184801543 M^(2//5) L^2 
 1.1300889873364839 
 0.038138404380536464 M^(9//2) L^(2//3) 
 1.1048123762940634 
 -0.6984190957741953 M^1 L^1 

julia> v .^ 2
5-element Vector{Quantity{Float64}}:
 0.6805558972463164 M^(4//5) L^4 
 1.2771011192991997 
 0.001454537888693323 M^9 L^(4//3) 
 1.220610386812535 
 0.48778923334204455 M^2 L^2
2 Likes
5

Okay this is probably my worst instance of Yak shaving yet, but I have now turned this into a new package:

Rather than using a Dict{Symbol,Bool} as above, it uses a fixed struct to store dimensions. There are 7 fields in this struct, corresponding to the 7 SI units.

Here’s the README:

This defines a simple statically-typed Quantity type for Julia.
Physical dimensions are stored as a value, as opposed to a parametric type, as in Unitful.jl.
This is done to allow for calculations where physical dimensions are not known at compile time.

Performance

DynamicUnits can greatly outperform Unitful
when the compiler cannot infer dimensions in a function:

julia> using BenchmarkTools, DynamicUnits; import Unitful

julia> dyn_uni = Quantity(0.2, mass=1, length=0.5, amount=3)
0.2 𝐋 ¹ᐟ² 𝐌 ¹ 𝐍 ³

julia> unitful = convert(Unitful.Quantity, dyn_uni)
0.2 kg m¹ᐟ² mol³

julia> f(x) = x ^ rand(1:10) * 0.3;

julia> @btime f($dyn_uni);
  80.449 ns (0 allocations: 0 bytes)

julia> @btime f($unitful);
  29.666 μs (42 allocations: 1.91 KiB)

(Note the μ and n.)
Here, the DynamicUnits quantity object allows the compiler to build a function that is type stable,
while the Unitful quantity object, which stores its dimensions in the type, requires type inference at runtime.

However, if the dimensions in your function can be inferred by the compiler,
then you are better off using Unitful:

julia> g(x) = x ^ 2 * 0.3;

julia> @btime g($dyn_uni);
  56.317 ns (0 allocations: 0 bytes)

julia> @btime g($unitful);
  1.958 ns (0 allocations: 0 bytes)

While both of these are type stable,
because Unitful parametrizes the type on the dimensions, functions can specialize
to units and the compiler can optimize away units from the code.

Usage

You can create a Quantity object with a value and keyword arguments for the powers of the physical dimensions
(mass, length, time, current, temperature, luminosity, amount):

julia> x = Quantity(0.3, mass=1, length=0.5)
0.3 𝐋 ¹ᐟ² 𝐌 ¹

julia> y = Quantity(10.2, mass=2, time=-2)
10.2 𝐌 ² 𝐓 ⁻²

Elementary calculations with +, -, *, /, ^, sqrt, cbrt are supported:

julia> x * y
3.0599999999999996 𝐋 ¹ᐟ² 𝐌 ³ 𝐓 ⁻²

julia> x / y
0.029411764705882353 𝐋 ¹ᐟ² 𝐌 ⁻¹ 𝐓 ²

julia> x ^ 3
0.027 𝐋 ³ᐟ² 𝐌 ³

julia> x ^ -1
3.3333333333333335 𝐋 ⁻¹ᐟ² 𝐌 ⁻¹

julia> sqrt(x)
0.5477225575051661 𝐋 ¹ᐟ⁴ 𝐌 ¹ᐟ²

julia> x ^ 1.5
0.1643167672515498 𝐋 ³ᐟ⁴ 𝐌 ³ᐟ²

Each of these values has the same type, thus obviating the need for type inference at runtime.

Furthermore, we can do dimensional analysis automatically:

julia> x + 3 * x
1.2 𝐋 ¹ᐟ² 𝐌 ¹

julia> x + y
INVALID

We can see the second one has valid(quantity) == false. This doesn’t throw an error by default, as it allows for stable return values.

The dimensions of a Quantity can be accessed either with dimension(quantity) for the entire Dimensions object:

julia> dimension(x)
𝐋 ¹ᐟ² 𝐌 ¹

or with umass, ulength, etc., for the various dimensions:

julia> umass(x)
1//1

julia> ulength(x)
1//2

Finally, you can strip units with ustrip:

julia> ustrip(x)
0.2

Units

Despite the name, DynamicUnits does not actually work with units. Instead, it works with dimensions.
You can use Unitful to parse units, and use the DynamicUnits->Unitful extension for conversion:

julia> using Unitful: Unitful, @u_str

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

julia> y = convert(DynamicUnits.Quantity, x)
500.0 𝐋 ¹ 𝐓 ⁻¹

julia> y2 = y^2 * 0.3
75000.0 𝐋 ² 𝐓 ⁻²

julia> x2 = convert(Unitful.Quantity, y2)
75000.0 m² s⁻²

julia> x^2*0.3 == x2
true

Vectors

There is not a separate class for vectors, but you can create units
like so:

julia> randn(5) .* Dimensions(mass=2/5, length=2)
5-element Vector{Quantity{Float64}}:
 -0.6450221578668845 𝐋 ² 𝐌 ²ᐟ⁵
 0.4024829670050946 𝐋 ² 𝐌 ²ᐟ⁵
 0.21478863605789672 𝐋 ² 𝐌 ²ᐟ⁵
 0.0719774550969669 𝐋 ² 𝐌 ²ᐟ⁵
 -1.4231241943420674 𝐋 ² 𝐌 ²ᐟ⁵

Because it is type stable, you can have mixed units in a vector too:

julia> v = [Quantity(randn(), mass=rand(0:5), length=rand(0:5)) for _=1:5]
5-element Vector{Quantity{Float64}}:
 2.2054411324716865 𝐌 ³
 -0.01603602425887379 𝐋 ⁴ 𝐌 ³
 1.4388184352393647 
 2.382303019892503 𝐋 ² 𝐌 ¹
 0.6071392594021706 𝐋 ⁴ 𝐌 ⁴
15 Likes
6

Thanks for doing this! I plan to test it out and recommend it over unitful for anything in SciML. I’ve been waiting for something like this for a long time and never got around to it, and so I’m very happy to see it finally exist.

3 Likes
7

Awesome, let me know if you have any suggestions for it :slightly_smiling_face:

1 Like
8

Just added a Unitful.jl extension too. So one can use Unitful.jl as a user-facing interface with its catalogue of units, and then DynamicUnits.jl for internal calculations:

julia> using Unitful: Unitful, @u_str

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

julia> y = convert(DynamicUnits.Quantity, x)
500.0 𝐋^1 𝐓^(-1)

julia> y2 = y^2 * 0.3
75000.0 𝐋^2 𝐓^(-2)

julia> x2 = convert(Unitful.Quantity, y2)
75000.0 m² s⁻²

julia> x^2*0.3 == x2
true

(It uses Unitful.upreferred to reduce a quantity to a standard form before stripping units, which seems to be fairly robust)

3 Likes
9

Awesome package, thanks! As long as it is not registered, we may wonder for a better meaningful name, eg I think of DimensionalUnits as you mention it works with Dimensions?

2 Likes
10

Thanks!

I am open to name suggestions.

To give backstory, here’s my original motivation for “DynamicUnits”:

Dynamic:

  1. Unitful.jl has dimensions as well; what differentiates this package is that the dimensions/units can change during runtime. “Dynamic” seemed like a good way to describe this.
  2. DynamicExpressions.jl is the cousin package and basically does the same thing for equations: it creates type stable symbolic expressions that can change during runtime without additional evaluations.

Units:

  1. Though it deals with physical dimensions, so long as you use a standardized unit system (such as I use when converting to/from Unitful - with upreferred), it is essentially the same thing. i.e., I guess I would argue that a “unit” is a physical dimension simply attached to a normalization factor.
    • It should be noted if not obvious from above that DynamicUnits tracks the quantity value as well as the dimensions. So you can totally convert to a unit system afterwards without losing any information.
  2. Subjectively, “dimension” seems like overloaded of a term, whereas “Unit” is associated to physical quantities. It also perhaps makes it easier “DynamicUnits” join into the Unitful.jl ecosystem of packages (which are pretty much all *Unit*.jl)

What do you think?

11

Regarding the name, I was thinking DimensionalQuantities.jl makes more sense for what it currently offers, given that there are no actual units. Obviously, this could evolve into something bigger though.

Thanks for putting this together.

2 Likes
12

That’s a good idea too. Actually I might change it to that…

I guess if the name is “Units,” people might assume it stores a catalog of units (maybe units that associated to dynamical systems…), like other packages in the Unitful.jl ecosystem. Whereas DynamicQuantities.jl is more obvious that it is an entire quantity type altogether.

Will think more about this

13

Given that a quantity is a more fundament concept for which units are based, it makes sense to me to have a package that defines a quantity type(s) that a separate unit package would build upon.

The use of “units” in the title is a double edge sword. It would make it easy for people to find, but I think most users would expect to find something that propagates units like Unitful or Pint in python

1 Like
14

Good point!

15

Actually I think I misread your suggestion initially. I had read it as “DynamicQuantities.jl” (which I do quite like and think I might switch to…). “DimensionalQuantities.jl” sounds like a bit of a tautology, maybe? But I think the “Quantities” part of the title is good.

16

How do error messages look like with your package? Less monstrous than those produced by using Unitful ?

17

I don’t do much for errors. The one thing I’d note is that instead of Unitful.DimensionError, I return the quantity as normal and set the (q::Quantity).valid flag to false. This is to remain consistent with the theme of type stability. That’s it though; it’s a pretty lightweight package: only ~300 SLOC.

(MethodError will still occur if a method is not defined for the Quantity type.)

This error return inspired by GitHub - JuliaPreludes/Try.jl: Zero-overhead and debuggable error handling. From the discussion on Performance of hasmethod vs try-catch on MethodError it seems like try-catch on a DimensionError is bound to be slower than tracking an extra value for errors in calculation.

18

This looks great! One question I have is whether Dimensions should have a type parameter which would let users with simple units use an Int8 per dimension rather than a SimpleRatio{Int}. Doing this would shrink a Dimensions object down from 896 bits to 56 bits (64 with padding), and would also let you compare Dimensions in a single instruction. For similar reasons, I think that Dimmensions should probably default to SimpleRatio{Int8}. I don’t think people commonly need more than 127 for their SI units

1 Like
19

If there is interest I would definitely welcome a PR on this. It may make the code a bit heavier but I think it is worth it:

I think exploring other types for the dimensions struct is a great idea in general. Even switching to @tim.holy’s SimpleRatio{Int} results in massive speedups compared to Rational{Int}. Now a benchmark against a compiled Unitful.jl method is only a 3x factor slower:

julia> dyn_uni = Quantity(0.2, mass=1, length=0.5, amount=3)
0.2 𝐋 ¹ᐟ² 𝐌 ¹ 𝐍 ³

julia> unitful = convert(Unitful.Quantity, dyn_uni)
0.2 kg m¹ᐟ² mol³

julia> g(x) = x ^ 2 * 0.3;

julia> @btime g($dyn_uni);
  6.083 ns (0 allocations: 0 bytes)

julia> @btime g($unitful);
  1.958 ns (0 allocations: 0 bytes)

And the non-compilable version is now >3000x faster than Unitful:

julia> f(x, i) = x ^ i * 0.3;

julia> @btime f($dyn_uni, i) setup=(i=rand(1:10));
  9.384 ns (0 allocations: 0 bytes)

julia> @btime f($unitful, i) setup=(i=rand(1:10));
  29.667 μs (42 allocations: 1.91 KiB)

(Note the μ and n.)

20

The package looks really nice - great initiative. I personally like the name DimensionalQuantities.jl, I don’t think it’s much of a tautology.

1 Like