SymbolicML / DynamicQuantities.jl

Efficient and type-stable physical quantities in Julia
https://symbolicml.org/DynamicQuantities.jl/dev/
Apache License 2.0
132 stars 17 forks source link

Create `AutoFloat` type for units #66

Open MilesCranmer opened 11 months ago

MilesCranmer commented 11 months ago

Attempt 2 at #30. cc @gaurav-arya - let me know what you think. This should fixes #65.

The following behavior occurs:

Input Resulting Numeric Type
Float16(1.0)u"km/s" Float16 type
1f0u"km/s" Float32 type
1.0u"km/s" Float64 type
big(1.0)u"km/s" BigFloat type

I also do:

Input Resulting Numeric Type
1u"km/s" Float64 type

The reason I convert integers to float is because of type instability. I don't want some units to convert to float for the sake of precision issues (like femtometers) while others can stay in integral form. So for integers the user will have to be explicit.

TODO remaining:

github-actions[bot] commented 11 months ago

Benchmark Results

main 3bc44636a87a4f... t[main]/t[3bc44636a87a4f...]
Quantity/creation/Quantity(x) 3.1 ± 0.1 ns 3.4 ± 0.099 ns 0.912
Quantity/creation/Quantity(x, length=y) 3.4 ± 0 ns 3.4 ± 0 ns 1
Quantity/with_numbers/*real 3.4 ± 0 ns 3 ± 0.1 ns 1.13
Quantity/with_numbers/^int 11.4 ± 3.6 ns 11.8 ± 4.3 ns 0.966
Quantity/with_numbers/^int * real 11.8 ± 4 ns 12.1 ± 4.3 ns 0.975
Quantity/with_quantity/+y 7.7 ± 0.1 ns 7.1 ± 0 ns 1.08
Quantity/with_quantity//y 4 ± 0 ns 3.4 ± 0 ns 1.18
Quantity/with_self/dimension 1.7 ± 0 ns 1.7 ± 0 ns 1
Quantity/with_self/inv 3.6 ± 0.3 ns 3.4 ± 0.2 ns 1.06
Quantity/with_self/ustrip 1.7 ± 0 ns 1.7 ± 0 ns 1
QuantityArray/broadcasting/multi_array_of_quantities 0.225 ± 0.019 ms 0.227 ± 0.027 ms 0.993
QuantityArray/broadcasting/multi_normal_array 0.0818 ± 0.0049 ms 0.0758 ± 0.0011 ms 1.08
QuantityArray/broadcasting/multi_quantity_array 0.255 ± 0.0021 ms 0.255 ± 0.0015 ms 1
QuantityArray/broadcasting/x^2_array_of_quantities 0.0464 ± 0.0059 ms 0.0444 ± 0.0033 ms 1.05
QuantityArray/broadcasting/x^2_normal_array 9.3 ± 1.5 μs 9.4 ± 1.7 μs 0.989
QuantityArray/broadcasting/x^2_quantity_array 10.3 ± 0.9 μs 10.4 ± 1.1 μs 0.99
QuantityArray/broadcasting/x^4_array_of_quantities 0.136 ± 0.0073 ms 0.137 ± 0.0066 ms 0.995
QuantityArray/broadcasting/x^4_normal_array 0.069 ± 0.0009 ms 0.0719 ± 0.0032 ms 0.96
QuantityArray/broadcasting/x^4_quantity_array 0.098 ± 0.011 ms 0.106 ± 0.0026 ms 0.922
time_to_load 0.17 ± 0.0004 s 0.178 ± 0.0008 s 0.957

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).

devmotion commented 11 months ago

This seems too complicated to me. Isn't it possible to just convert apply float to parsed quantities as suggested in #65 but otherwise not modify the user input?

MilesCranmer commented 11 months ago

This seems too complicated to me. Isn't it possible to just convert apply float to parsed quantities as suggested in #65 but otherwise not modify the user input?

I am not sure how to do it any other way. Do you have an implementation in mind?

devmotion commented 11 months ago

To me it seems the main problem is that the current implementation of the u_str macro always returns a Quantity. If no values are given, I think it should just return a Dimensions. This would also match the behaviour of Unitful:

julia> typeof(Unitful.u"m")
Unitful.FreeUnits{(m,), 𝐋, nothing}

julia> typeof(Unitful.u"1f0m")
Unitful.Quantity{Float32, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}

julia> typeof(1f0Unitful.u"m")
Unitful.Quantity{Float32, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}

julia> typeof(Unitful.u"1f0")
Float32

I think this would be more intuitive than the current behaviour of DynamicQuantities:

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

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

julia> typeof(1f0DynamicQuantities.u"m")
DynamicQuantities.Quantity{Float64, DynamicQuantities.Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}

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

Edit: I retract this comment partly. I guess these inconsistencies with Unitful are unavoidable since DynamicQuantities does not store the unit in the type domain. However, I don't think DynamicQuantities should try to convert everything to the default type.

MilesCranmer commented 11 months ago

These conversions are to make the return value of uparse type stable. In Unitful you can even do u"(m, km, cm)". However in DynamicQuantities this is not allowed since it would make uparse return ::Any.

I think it should just return a Dimensions

The operator * does not exist between Number and Dimensions which is on purpose. This is because it would necessitate choosing a default AbstractQuantity.

MilesCranmer commented 11 months ago

Also keep in mind (I see your edit mentioned this) that all units and constants have numerical values in SI base units, so they are stored as Quantity rather than just Dimensions.

devmotion commented 11 months ago

These conversions are to make the return value of uparse type stable

It's a macro, so the string is replaced with whatever expression you want to replace it with. And then the resulting function/code/etc. is compiled by the Julia compiler. I don't see how the conversions help type stability.

As a quick demonstration I replaced the uparse function with just return eval(Meta.parse(s)) (BTW I wonder if one could omit eval and just replace the m etc. Symbols with $(Units.m) etc., such that all evaluation is only performed after the macro has been applied). There are no type instability issues:

julia> using DynamicQuantities, Test

julia> julia> f(x) = x * u"3m/s"
f (generic function with 1 method)

julia> @inferred f(3)
9.0 m s⁻¹

julia> @inferred f(1.0)
3.0 m s⁻¹

julia> @code_warntype f(2)
MethodInstance for f(::Int64)
  from f(x) @ Main REPL[18]:1
Arguments
  #self#::Core.Const(f)
  x::Int64
Body::Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}
1 ─ %1 = (x * 3.0 m s⁻¹)::Core.PartialStruct(Quantity{Float64, Dimensions{DynamicQuantities.FixedRational{Int32, 25200}}}, Any[Float64, Core.Const(m s⁻¹)])
└──      return %1

julia> @code_lowered f(1.0)
CodeInfo(
1 ─ %1 = x * 3.0 m s⁻¹
└──      return %1
)

This is because it would necessitate choosing a default AbstractQuantity.

To be precise, it would only require to choose a default quantity type for Dimensions. And in that case it might be fine to choose Quantity since it is the default quantity type anyway and is already based on Dimensions?

MilesCranmer commented 11 months ago

I guess using @u_str is type stable even if uparse is type unstable. Hmm...

However I do like having u"1" being a Quantity. It returns the unit Quantity type with default dimensions and types, so you don’t have to write it out manually. If this returns just a literal 1, and you would call it at compile time, what is the point of it?

We need AutoFloat regardless because any quantities (i.e., constants or units) store numerical values, unlike Unitful, where you can have separate symbolic dimensions in the type. So we need a way to demote to Float16.

devmotion commented 11 months ago

I still think that it would be preferable to change the macro and parsing in more fundamental ways. But in any case, I don't think adding a new numeric type is a good idea (I've seen too many things go wrong in different packages). The much simpler approach would be to make the "units" quantities with Boolean values. Bools are frequently used to avoid undesired promotions since they are the least invasive default number types in Julia:

julia> 1.0 * true
1.0

julia> 2 * true
2

julia> Float16(3) * true
Float16(3.0)

julia> 4f0 * true
4.0f0
MilesCranmer commented 11 months ago

Unlike Bool I think we want to convert integers to floats, otherwise the user could write 1u"cm" and get hit with a rounding error which seems unpleasant.

devmotion commented 11 months ago

I'm not sure if that's a good idea. To me it seems rather surprising to not respect the types chosen by a user (which also lead to #65 and is inconsistent with Unitful).

If you really want to prevent this case, you could specialize *(::Int, ::Quantity{Bool}) etc. (I'm a bit worried that specializing *(::Integer, ::Quantity{Bool}) creates invalidations but this should be checked, of course).

MilesCranmer commented 11 months ago

The user could still convert to a Quantity{Int} afterwards. It’s not so much it doesn’t respect their choice, it’s more: the units and constants are already floating point numerical values, and promote accordingly. Writing 1u"cm" is multiplying an integer with a float.

MilesCranmer commented 11 months ago

Btw with AutoFloat, you do get the behavior of 1f0"cm" being a Float32. (See the table above)

gaurav-arya commented 11 months ago

To clarify the key problem: the question is, when the user writes something like u"fm", what should the resulting quantity look like? Since the quantity will be represented in SI, a 1e-15 has to be written somewhere, and thus a floating point type needs to be chosen. After understanding this, I see why the current u-str behaviour is what it is. In theory users could always use the proper constructor Quantity to specify the value type (I believe? Or does that also promote to Float64?) However, the u-str API is the most convenient way to create a quantity right now, and furthermore the promotion behaviour would be surprising to those who do not understand the internals of DynamicQuantities.jl. This raises the question: can we make the recommended / easiest to use API for constructing quantities more generic with respect to the value type?

It's a tricky problem and I don't have a perfect solution, but hopefully the thoughts below are helpful.

  1. The simplest solution from an implementation perspective would be for the recommend API to require users to specify both the value and the units of the quantity at once. This would not require AutoFloat. The drawback is that the u-str approach exactly as is would no longer be OK, we'd have to tweak the API. One API here (which takes inspiration from option 2.) would be to enforce expressions of the form u"x * m" in u-strs and syntactically extract x and use it as the value. (Edit: instead of doing it syntactically, a really cute way of enforcing u"x * m" would be to use the approach described in 2 below, but throw an error if a u-str parses to a quantity that is still backed by an AutoFloat.)

  2. Alternatively, one could view the form x * u"m" as essentially a fancy constructor for quantities, made possible with the AutoFloat logic. In this solution (the one currently implemented by this PR), I would recommend strongly discouraging creating plain unit strings u"m" (even if x happens to be 1.0), as the form x * u"m" is how units should be constructed (or even u"x * m", see the last sentence of this paragraph). In this perspective, I think AutoFloat has a nice interpretation: it is fancy magic for getting this Unitful-like constructor syntax to work with this package, but if quantity are constructed as recommended then it should dissapear immediately for subsequent calculations. Furthermore, if the user writes u"x * m" with x inside the u-str, no AutoFloats would appear at runtime whatsoever, only in the construction at compile time.

Edit: lastly, regarding @devmotion's concern about adding a new numerical type, I very much share it: it is indeed a somewhat tricky problem. Essentially, it boils down to there not being an AbstractWrapperNumber that you can just subtype and overload a few methods of OOP-style (x-ref https://github.com/SimonDanisch/AbstractNumbers.jl).

However, two points regarding AutoFloat: 1) This package has to solve this problem anyway with Quantity. With the right design (nontrivial, admittedlly...), we would be able to write overloads that were figured out for Quantity with any other numerical types too. And more importantly, 2) as explained above, we would want to absolutely minimize (i.e. relegate only to construction logic) or even eliminate the use of AutoFloat in runtime: if it exists it ought to be a tool for construction, and not go any further than that. It is true than in option 2. we may not be able to 100% prevent users from working with a plain AutoFloat-backed unit-string like u"m" even if it is strongly discouraged. But perhaps error hints on MethodErrors with AutoFloat's explaining this common user error could make this even more robust.

Edit 2: in fact, as a middle ground, what if we intentionally only defined * on AutoFloat, and threw method errors with error hints for everything else?

MilesCranmer commented 11 months ago

Thanks for sharing all of your thoughts on this!

The one comment I will make is that I think we should try to avoid an API that requires things like u"x * m". I think over-reliance on macros is not a great code pattern and introduces brittleness and new bugs (e.g., what if a user wants to have x as a variable?). Also we would still want to allow the user to do things like:

using DynamicQuantities.Units: m, kg
using DynamicQuantities.Constants: c

q = 0.3f0 * kg/m^3 * c^2

which is just normal Julia code (which, btw, would create a Quantity{Float32} after this PR). Those m, kg, and c need to have some type after all, so the solution should be compatible with this as well.

Keep in mind this AutoFloat proposal is just a thin wrapper around Float64 so I wouldn't expect it to have any runtime penalties as the compiler would inline stuff. But yeah I 100% agree that if we go this route we should try to force the user to never use AutoFloat in any computations; it's just there for storing the physical unit/constant's value – but should get converted immediately. The idea of throwing easy-to-understand errors for any other pattern is a good idea.

MilesCranmer commented 10 months ago

Just to clarify something – the promotion rules here turn out to be the same as used in irrational numbers. e.g.,

julia> promote(π, 1) |> typeof
Tuple{Float64, Float64}

julia> promote(π, 1f0) |> typeof
Tuple{Float32, Float32}