DynamicQuantities.jl v0.7.0: efficient and type-stable physical quantities

Happy to share DynamicQuantities.jl v0.7.0!

image

Dev Build Status Coverage Aqua QA

The package has had some major changes over the last few versions, so I think it deserved a new post.

DynamicQuantities defines a simple statically-typed Quantity type for storing physical units.

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.

Changes since the last post (described more below)

  • Performance improvements (via better compiler inlining)
  • Units, and unit parsing via the @u_str macro
  • Physical constants, and constant parsing via the Constants.* prefix (also in @u_str)
  • SymbolicDimensions for working with symbolic units and constants
    • Standard Dimensions will eagerly convert to SI units; this avoids it.
    • Symbolic unit/constant parsing via the @us_str macro.
  • QuantityArray <: AbstractArray for efficiently storing arrays of quantities that have the same units
  • Extensions for ScientificTypes.jl, LinearAlgebra.jl, as well as an extension to convert to/from Unitful.jl quantities
  • AbstractQuantity and AbstractDimensions for extending behavior or defining custom spaces of physical dimensions

Thanks to contributions from Gaurav Arya @odow @Oscar_Smith @j-fu, and suggestions/feedback from @non-Jedi @ChrisRackauckas @mcabbott and many others in the last thread and on GitHub!

Performance

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

julia> using BenchmarkTools, DynamicQuantities; import Unitful

julia> dyn_uni = 0.2u"m^0.5 * kg * mol^3"
0.2 mĀ¹įŸĀ² kg molĀ³

julia> unitful = convert(Unitful.Quantity, dyn_uni)
0.2 kg mĀ¹įŸĀ² molĀ³

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

julia> @btime f($dyn_uni, 1);
  2.750 ns (0 allocations: 0 bytes)

julia> @btime f($unitful, 1);
  2.718 Ī¼s (30 allocations: 1.34 KiB)

Note the Ī¼ and n: this is a 1000x speedup. Here, the DynamicQuantities 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 can get better speeds with Unitful:

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

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

julia> @btime g($unitful);
  1.500 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.

Aside: The compiler seems to be pretty good at inlining things especially with the recent update (thanks to Gaurav Arya for helping get this through!), so this performance gap seems to have shrunk for even statically typed calculations. This above calculation used to be 5x in favor of Unitful. However, this depends on compiler constant propagation so it is calculation dependent how big this gap would be.

Types

At the heart of the package is just two immutable structs:

struct Dimensions{R<:Real} <: AbstractDimensions{R}
    length::R
    mass::R
    time::R
    current::R
    temperature::R
    luminosity::R
    amount::R
end

struct Quantity{T,D<:AbstractDimensions} <: AbstractQuantity{T,D}
    value::T
    dimensions::D
end

The R type here is typically a rational-like number (the default is an internal FixedRational type ā€“ a rational number with fixed denominator which gives much faster operations than Rational).

Whatā€™s nice about the abstract interface is you can just write a custom AbstractDimensions type with the physical dimensions you want to use, and via the use of @oxinaboxā€™s Tricks.jl, static_fieldnames is used to compile all of the unit propagation methods at first call. So e.g., struct MyDimensions{R} <: AbstractDimensions{R}; length::R; time::R; end would work out of the box! (Quantity(0.3, MyDimensions, length=1, mass=-1) would then be 0.3 m/s, and you could calculate away)

Generic usage

You can create a Quantity object by using the convenience macro u"...":

julia> x = 0.3u"km/s"
300.0 m sā»Ā¹

julia> room_temp = 100u"kPa"
100000.0 mā»Ā¹ kg sā»Ā²

This supports a wide range of SI base and derived units, with common prefixes.

You can also construct values explicitly with the Quantity type, with a value and keyword arguments for the powers of the physical dimensions (mass, length, time, current, temperature, luminosity, amount):

julia> x = Quantity(300.0, length=1, time=-1)
300.0 m sā»Ā¹

Elementary calculations with +, -, *, /, sqrt, cbrt, abs are supported, and ^ will use rationalize to get a reasonable power from exponentiation:

julia> x * y
12600.0 m kg sā»Ā¹

julia> x / y
7.142857142857143 m kgā»Ā¹ sā»Ā¹

julia> x ^ 3
2.7e7 mĀ³ sā»Ā³

julia> x ^ -1
0.0033333333333333335 mā»Ā¹ s

julia> sqrt(x)
17.320508075688775 mĀ¹įŸĀ² sā»Ā¹įŸĀ²

julia> x ^ 1.5
5196.152422706632 mĀ³įŸĀ² sā»Ā³įŸĀ²

Each of these values has the same type, which means we donā€™t need to perform type inference at runtime.

Furthermore, we can do dimensional analysis by detecting DimensionError:

julia> x + 3 * x
1.2 mĀ¹įŸĀ² kg

julia> x + y
ERROR: DimensionError: 0.3 mĀ¹įŸĀ² kg and 10.2 kgĀ² sā»Ā² have incompatible dimensions

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

julia> dimension(x)
mĀ¹įŸĀ² kg

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

Constants

There are a variety of physical constants accessible
via the Constants submodule:

julia> Constants.c
2.99792458e8 m sā»Ā¹

These can also be used inside the u"..." macro:

julia> u"Constants.c * Hz"
2.99792458e8 m sā»Ā²

For the full list, see the docs.

Symbolic Units

You can also choose to not eagerly convert to SI base units, instead leaving the units as the user had written them. For example:

julia> q = 100us"cm * kPa"
100.0 cm kPa

julia> q^2
10000.0 cmĀ² kPaĀ²

You can convert to regular SI base units with expand_units:

julia> expand_units(q^2)
1.0e6 kgĀ² sā»ā“

This also works with constants:

julia> x = us"Constants.c * Hz"
1.0 Hz c

julia> x^2
1.0 HzĀ² cĀ²

julia> expand_units(x^2)
8.987551787368176e16 mĀ² sā»ā“

This dimensions type works a bit differently as it stores all the dimensions in a sparse vector (source for the curious). All unit calculations are performed as operations on this sparse vector (this is convenient for user interfaces, and still reasonably fast, but not as fast as the regular immutable Dimensions).

Arrays

For working with an array of quantities that have the same dimensions, you can use a QuantityArray:

julia> ar = QuantityArray(rand(3), u"m/s")
3-element QuantityArray(::Vector{Float64}, ::Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}):
 0.2729202669351497 m sā»Ā¹
 0.992546340360901 m sā»Ā¹
 0.16863543422972482 m sā»Ā¹

This QuantityArray is a subtype <:AbstractArray{Quantity{Float64,Dimensions{...}},1},
meaning that indexing a specific element will return a Quantity:

julia> ar[2]
0.992546340360901 m sā»Ā¹

julia> ar[2] *= 2
1.985092680721802 m sā»Ā¹

julia> ar[2] += 0.5u"m/s"
2.485092680721802 m sā»Ā¹

Which performs dimension checks.

This has a custom broadcasting interface which allows the compiler to avoid redundant dimension calculations and dimensional analysis, relative to if you had simply used an array of quantities:

julia> f(v) = v^2 * 1.5;

julia> @btime $f.(xa) setup=(xa = randn(100000) .* u"km/s");
  109.500 Ī¼s (2 allocations: 3.81 MiB)

julia> @btime $f.(qa) setup=(xa = randn(100000) .* u"km/s"; qa = QuantityArray(xa));
  50.917 Ī¼s (3 allocations: 781.34 KiB)

So we can see the QuantityArray version saves on both time and memory.

Unitful

DynamicQuantities allows you to convert back and forth from Unitful.jl:

julia> import Unitful; import DynamicQuantities

julia> x = 0.5Unitful.u"km/s"
0.5 km sā»Ā¹

julia> y = convert(DynamicQuantities.Quantity, x)
500.0 m sā»Ā¹

julia> y2 = y^2 * 0.3
75000.0 mĀ² sā»Ā²

julia> x2 = convert(Unitful.Quantity, y2)
75000.0 mĀ² sā»Ā²

julia> x^2*0.3 == x2
true

Other

One other thing to mention, which is actually part of the reason I started working on this package: SymbolicRegression.jl now supports dimensional constraints as part of equation discovery searches! This feature required the creation of this package, as expressions are runtime generated so Unitful.jl could not be used efficiently. Example here: Examples Ā· SymbolicRegression.jl

42 Likes

Thanks for the package! Itā€™ll also be useful outside of the diffeq ecosystem sometimes.

I wonder how compatible its interface is with Unitful? That is, can one take code, replace using Unitful with DynamicQuantities, and be done? Or, is this among future goals?

Also, just a minor comment:
The name clearly conveys what this package does, but the short description doesnā€™t really differentiate it from Unitful.

fast and lightweight physical quantities

simple statically-typed type for storing physical units

Both are applicable to Unitful as-is. Maybe, highlight dynamic/runtime/untyped nature more?

1 Like

That is, can one take code, replace using Unitful with DynamicQuantities, and be done?

I think this is mostly there? The same core symbols: @u_str, uparse, and ustrip can be used the same way (mostly; I donā€™t allow u"(m, s)" for example as it would require type instability), and @u_str creates a Quantity. This Quantity type in both packages forwards all the basic mathematical methods you would expect. You can also dispatch on Quantity{Float64} in both packages to expect quantities with a Float64 value type. But please raise an issue if there are any major incompatibilities Iā€™m missing! I would like it to be compatible aside from the major design differences ā€“

The main difference which we are stuck with is from how they treat dimensions at the type level (if that is what you are asking). So if you are dispatching on dimensions, like

julia> using Unitful

julia> x = 1.5u"m/s"
1.5 m sā»Ā¹

julia> typeof(x)
Quantity{Float64, š‹ š“ā»Ā¹, Unitful.FreeUnits{(m, sā»Ā¹), š‹ š“ā»Ā¹, nothing}}

and then dispatching on <:Quantity{Float64, š‹ š“ā»Ā¹}, then you would have to change that code to look at the values of the quantity instead. Because DynamicQuantities.jl does this:

julia> using DynamicQuantities

julia> x = 1.5u"m/s"
1.5 m sā»Ā¹

julia> typeof(x)
Quantity{Float64, Dimensions{FixedRational{Int32, 25200}}}

This is the same type regardless of the units used. The dimensions are stored in dimension(x) which is a simple immutable struct:

julia> Base.propertynames(dimension(x))
(:length, :mass, :time, :current, :temperature, :luminosity, :amount)

which are all rationals.

Good point, thanks! Any suggestions for wording it would be much appreciated. Maybe something like DynamicQuantities defines a simple type-consistent Quantity for storing physical units?

This is very nice indeed. I do suspect however that this is not for large-scale simulations? There is a penalty to pay when storing and manipulating a large number of numbers with units?

1 Like

Thanks. So If you are confident that there are no type instabilities due to failure to infer units during compilation (which can really hurt performance), or even things like multiple units in a single array (which IIRC might promote to Array{Any}?), then yes Unitful will be faster.

How much faster depends on the calculation and how much inlining the compiler does.

Sometimes the compiler can propagate unit values pretty well by itself, without needing to encode everything in the type (as seen in the above comparison). And the QuantityArray means it might choose to process the units only once for the entire array. For storage costs, this would be a single Dimensions for the whole array ā€” there is not a separate Dimensions stored for each number. So only 32 bytes extra per entire array.

My gut feeling is that if this is a highly tuned and type stable simulation with a few expensive kernels that are called repeatedly, then I might stick with Unitful as it would completely eliminate the cost. But since the API is pretty similar you might as well try and see :wink:

The other pro is you can take full advantage of precompilation, since you donā€™t need to compile the simulation for all possible units, only once for a generic Quantity.

3 Likes

Not inlining but branch prediction. The switch statements of ā€œif this unitā€ etc. are actually ā€œfreeā€ when running in a loop, because your CPU does speculative execution and so it just starts executing one branch under the assumption that itā€™s correct and if itā€™s correct thereā€™s no cost. When run in a loop, this has a process of effectively ā€œlearningā€ the branches to predict (internal on the CPU, some cool stuff) and so the branches predict better and the cost of having the branches is effectively zero. See:

Given this fact, I think way too many people overestimate the utility of putting this unit information into the type domain. That said, the existence of branches do cause LLVM to omit SIMD, and so if you did have a loop that SIMDā€™d very well you would pay the price of that being eliminated because the runtime checks would interfere with the assumptions of SIMD.

But that said, for a lot of ODE/PDE models, the main cost is in the LU-factorizations and if that just ignores units, this shouldnā€™t be too noticable of a hit to the f evaluations and so it shouldnā€™t actually show up as a major contributor in ā€œmostā€ profiles. You can definitely find some cases where this matters (very small cases getting very good SIMD), but I think people should really go to this technique first before trying Unitful in almost every case because Unitful is a lot more work for a performance optimization that only exists in a minority of cases (especially because many cases with Unitful will hit Array{Any} and actually be orders of magnitude slower due to dynamic dispatching).

Moral of the story, DynamicQuantities.jl is really good.

9 Likes

I think the issue is that the additional overhead due to units is often due to running something like unit1 * unit2 = unit3 repeatedly in a loop? This computation isnā€™t a branch but rather some arithmetic operations that are redundantly performed each iteration.

Using a QuantityArray we could hope for this to be lifted out of the loop since the compiler should in theory be able to prove that unit1 and unit2 are loop invariant, but this only seems to happen in limited cases dependent upon inlining heuristics.

(Spitballing here, but if pure branch computations are free, I wonder if we could optimize unit1 * unit2 to be a pure branch computation in common cases via a lookup table-esque approachā€¦)

The compiler doesnā€™t lift it out in this case, but branch prediction can make the cost trivial anyways because there is no cost to well-predicted branches, and a simple structure will predict well. Again, what you lose in performance is just SIMD.

2 Likes

Some more updates as of v0.10

New typing system

There are now multiple quantity types which act identically, but let you work with different subtypes:

  • Quantity <: AbstractQuantity <: Number
    • The default quantity returned by 0.5u"m/s", which is subtyped to Number.
  • GenericQuantity <: AbstractGenericQuantity <: Any
    • Do you need to put units on some non-numerical type? Then you want GenericQuantity
  • RealQuantity <: AbstractRealQuantity <: Real
    • Many packages throughout Julia require Real input. Although physical quantities are not technically real, it can still be useful for compatibility to pass them as if they were. RealQuantity lets you do that.
  • All three of these are of course compatible with QuantityArray, and can take any AbstractDimensions to store the dimensions.
julia> x = 0.5u"m/s"
0.5 m sā»Ā¹

julia> x_re = RealQuantity(x)::Real
0.5 m sā»Ā¹

julia> y_gen = GenericQuantity("hello", length=1, time=-1)
(hello) m sā»Ā¹

We also get automatic promotion:

julia> x = RealQuantity(0.5u"m/s")
0.5 m sā»Ā¹

julia> y = x * (1 + 2im)
(0.5 + 1.0im) m sā»Ā¹

julia> typeof(y)
Quantity{ComplexF64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}

the behavior of which can be overloaded for any custom hierarchy of quantity types.

Thanks to Guarav Arya for his help on this.

New methods

Many other numerical methods which require dimensionless numbers, such as exp or sin, can now take a quantity as input (required because DynamicQuantities.Quantity can be dimensionless - it doesnā€™t automatically convert to Float64 due to the potential for type instability).

The overloaded method for dimensionless functions will simply check if the quantity is dimensionless or not, and then evaluate:

julia> x = 0.5u"m"
0.5 m

julia> y = 10u"Constants.Mpc"
3.085677581491367e23 m

julia> exp(x/y)
1.0

This means we can also have ranges of quantities now (via Base.rem)

julia> for x in 0u"km":1.5u"km":10u"km"
           @show x^2
       end
x ^ 2 = 0.0 mĀ²
x ^ 2 = 2.25e6 mĀ²
x ^ 2 = 9.0e6 mĀ²
x ^ 2 = 2.025e7 mĀ²
x ^ 2 = 3.6e7 mĀ²
x ^ 2 = 5.625e7 mĀ²
x ^ 2 = 8.1e7 mĀ²

(Although this should probably use a QuantityArrayā€¦ PRs always appreciated )

Better perf

QuantityArray performance seems pretty solid now, and lets you wrap arrays of arbitrary data with quantities via GenericQuantity. Hereā€™s summation on a regular array compared to an array of physical quantities:

julia> struct Coords
           x::Float64
           y::Float64
       end

julia> Base.:+(a::Coords, b::Coords) = Coords(a.x+b.x, a.y+b.y)

julia> Base.:*(a::Coords, b::Number) = Coords(a.x*b, a.y*b)

julia> Base.:*(a::Number, b::Coords) = Coords(a*b.x, a*b.y)

julia> @btime(
    sum(array),
    setup=(N=1000; array=[Coords(rand(), rand()) for i=1:N]
)
  828.000 ns (0 allocations: 0 bytes)
Coords(509.5939555433635, 494.56896926222123)

julia> @btime(
    sum(array),
    setup=(N=1000; array=QuantityArray([GenericQuantity(Coords(rand(), rand()), dimension(u"m/s")) for i=1:N])
)
  843.750 ns (0 allocations: 0 bytes)
(Coords(504.7718424677944, 482.0449323023619)) m sā»Ā¹

Thanks to all the contributors thus far :smiley:

20 Likes

More updates as of v0.13.0!

Unit conversion operator

There is now a dedicated infix operator for converting units, thanks to @j-fu:

julia> 5e-9u"m" |> us"nm"
5.0 nm

which is equivalent to uconvert(us"nm", 5e-9u"m").

Custom units

In addition to custom dimension spaces which you can define by inheriting from AbstractDimensions:

julia> struct CookiesAndMilk{R} <: AbstractDimensions{R}
           cookies::R
           milk::R
       end

julia> cookie_rate = Quantity(0.9, CookiesAndMilk(cookies=1, milk=-1))
0.9 cookies milkā»Ā¹

julia> total_milk = Quantity(103, CookiesAndMilk(milk=1))
103 milk

julia> total_cookies = cookie_rate * total_milk
92.7 cookies

You can also easily define custom units in SI dimensions (i.e., built-in Dimensions), thanks to @ven-k:

julia> @register_unit OneFiveV 1.5u"V"

which lets you use this in normal macros:

julia> x = us"OneFiveV"
1.0 OneFiveV

julia> x * 10u"A" |> us"W"
15.0 W

julia> 3us"V" |> us"OneFiveV"
2.0 OneFiveV

Note that any custom units registered with @register_unit are localized to a particular module and do not leak into the global namespace.

Angular units

There are now angular units such as rad, arcmin, deg, and solid angle unit sr. Just note that these assume the SI definition that 1 rad = 1 = 1 sr, but you can always track them as symbolic units like:

julia> Īø_moon = 31us"arcmin"
31.0 arcmin

julia> d_moon = 384_400us"km"
384400.0 km

julia> r_moon = d_moon * Īø_moon / 2 |> uexpand
1.7331701248721024e6 m

julia> r_moon |> us"Constants.R_earth"
0.2717376844000725 R_earth

Symbolic unit exports

The primitive symbolic units are now stored using a single integer, and only create allocations when you start doing operations with them. This means that all symbolic units have now been made available for direct import:

julia> using DynamicQuantities: SymbolicUnits as U

julia> U.km / U.yr
1.0 km yrā»Ā¹

Thanks to @devmotion for help in implementing this.

Other changes

You can now precompile statements made with us"..." or u"..." as they no longer rely on eval but on an internal parser. Thanks to @YingboMa for tips.

There is now a proper isapprox for arrays of quantities and QuantityArray which allows atol with unit specification. Thanks to @mike.ingold for this.

9 Likes

Great seeing this gaining more and more features, becoming more and more complete Unitful alternative for those cases when runtime units are required!

However, a bit unfortunate that these error-prone defaults proliferate :(
(same as Unitful, but unlike eg astropy)
Allows for completely meaningless conversions:

julia> (u"W" / 30us"arcmin^2") |> us"W"
393936.76200140937 W

I donā€™t think thereā€™s really a ā€œcorrectā€ answer here. 1 rad = 1 is the SI standard: Radian - Wikipedia, so this statement of u"W"/u"rad" == u"W" is 100% correct according to the international system of units, even though it might not be suitable for a particular task.

In astropy their choice to violate SI means that you get things like:

>>> from astropy import units as u

>>> theta = 0.5*u.rad

>>> distance = 100*u.km

>>> (theta*distance).to(u.km)
...
astropy.units.core.UnitConversionError: 'km rad' and 'km' (length) are not convertible

Itā€™s really just a choice someone has to make. Astronomy, as you know, has countless examples of nomenclature that completely violate SI. (Sometimes itā€™s very convenient though!)

Iā€™d rather stick 100% to SI though (even though I myself have an astronomy background). There are just too many silent bugs that can arise from mixing standards.

4 Likes

Itā€™s much easier to tell the computer to ignore radians when they are present (eg, see astropy so-called ā€œequivalenciesā€) than to determine what angular units should be in the quantity when there are none.
So, the astropyā€™s approach makes both usecases (want angular units present or ignored) possible and convienient, while Unitfulā€™s only allows one.

Hard to imagine any silent bugs arising from more strict conversion rules, with more constraints. But bugs in the other direction are clearly possible (see my last post).

Wait, why canā€™t you do this with the existing symbolic units? Like us"arcmin". Only u"arcmin" will automatically convert to SI equivalencies.

This is with ā€œsymbolic unitsā€:

Itā€™s not:

julia> typeof(u"W" / 30us"arcmin^2")
Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}

You need to write us"W" / 30us"arcmin^2" to keep it in symbolic units. With u"W" it will promote to regular SI units.

julia> typeof(us"W" / 30us"arcmin^2")
Quantity{Float64, SymbolicDimensions{DynamicQuantities.FixedRational{Int32, 25200}}}

Same result anyways:

julia> (us"W" / 30us"arcmin^2") |> us"W"
393936.7620014093 W

Because you doing explicit conversion. If you write it out it stays there.

julia> x = us"W" / 30us"arcmin^2"
0.03333333333333333 W arcminā»Ā²

julia> y = x * 10us"arcmin"
0.3333333333333333 W arcminā»Ā¹

julia> y * 0.5us"arcmin"
0.16666666666666666 W

Does this help?

A main/major point of using unitful values instead of unitless (= raw numbers) is that only unitfully-meaningful conversions are allowed. Clearly, W/arcmin^2 being convertible to W doesnā€™t help at all :)

I had this bug quite a few times (with Unitful, but same behavior here): forgot to multiply by an angle, got correct units in the end, but completely meaningless results.

I mean, W/arcmin^2 is by definition equivalent dimension to W in SI. I guess here you just mean itā€™s not convenient for astronomy (which is very fair!).

I guess you could try to do something like this if you want?

julia> struct AngleDimensions{R} <: AbstractDimensions{R}
           length::R
           mass::R
           time::R
           current::R
           temperature::R
           luminosity::R
           amount::R
           angle::R
       end

julia> y = Quantity(2.0, AngleDimensions(mass=1, length=2, time=-3))
2.0 mĀ² kg sā»Ā³

julia> x = Quantity(0.5, AngleDimensions(angle=2))
0.5 angleĀ²

julia> y / x
4.0 mĀ² kg sā»Ā³ angleā»Ā²

It will need converters from regular Dimensions but the operations themselves should work out of the box.

1 Like