PainterQubits / Unitful.jl

Physical quantities with arbitrary units
Other
594 stars 111 forks source link

This package assumes that `*` is commutative #607

Closed sostock closed 1 year ago

sostock commented 1 year ago

Methods like the following are generally wrong (consider multiplying a Quaternion by a Quantity{Quaternion}): https://github.com/PainterQubits/Unitful.jl/blob/e2aa6f8bb66ddf37adc0d33c7a40958dfcaaadaa/src/quantities.jl#L32

cstjean commented 1 year ago

Out of curiosity, is a unitful quaternion a meaningful quantity? What would be the units?

sostock commented 1 year ago

There is a package which depends on both Quaternions.jl and Unitful.jl and (according to its documentation) explicitly supports Unitful: QuaternionIntegrator.jl. In my opinion, this makes our treatment of multiplication downright dangerous.

Calling QuaternionIntegrator.rotate(q, v) with a unitless quaternion q and a unitful vector v leads to wrong results, because it calls the method discussed above:

julia> using QuaternionIntegrator, Quaternions, Unitful

julia> q = Quaternion(1,2,3,4);

julia> v = [1,2,3];

julia> rotate(q, v)
3-element Vector{Float64}:
 1.7999999999999998
 2.0
 2.5999999999999996

julia> rotate(q, v*u"m")
3-element Vector{Quantity{Float64, š‹, Unitful.FreeUnits{(m,), š‹, nothing}}}:
 1.0 m
 2.0 m
 3.0 m

I would assume that this also affects the actual usecases of QuaternionIntegrator. Unfortunately, the QuaternionIntegrator tests for unitful vectors only check that the output has the correct units: https://github.com/dronir/QuaternionIntegrator.jl/blob/2db6a2634e6726920d190dad4fcef4f2b4c184e6/test/runtests.jl#L48-L60.

The results of the test with and without units are not the same:

āˆ†t = 0.001 * u"s"
I = diagm([1.0, 1.0, 1.0]) * u"kg * m^2"
torque(q) = [0.1, 0.0, 0.0] * u"N * m"
q0 = Quaternion(1.0, 0.0, 0.0, 0.0)
w0 = [0.0, 0.0, 1.0] * u"1/s"

q1, w1 = integrate(q0, w0, I, āˆ†t, torque)

# result:
q1 = QuaternionF64(0.9999998750000023, 2.499999895833335e-8, -8.271805780871692e-28, 0.000499999979166667, false)
w1 = [9.999999999999998e-5 s^-1, 3.308722450212111e-24 s^-1, 0.9999999999999999 s^-1]
āˆ†t = 0.001
I = diagm([1.0, 1.0, 1.0])
torque(q) = [0.1, 0.0, 0.0]
q0 = Quaternion(1.0, 0.0, 0.0, 0.0)
w0 = [0.0, 0.0, 1.0]

q1, w1 = integrate(q0, w0, I, āˆ†t, torque)

# result:
q1 = QuaternionF64(0.9999998750000023, 2.4999997395833446e-8, 6.249999479166684e-12, 0.000499999979166667, false)
w1 = [0.00010000002499999401, 1.5624998339452638e-14, 0.9999999999999976]
brainandforce commented 1 year ago

Out of curiosity, is a unitful quaternion a meaningful quantity? What would be the units?

There's a nice way of interpreting quaternions in the framework of geometric algebra - the geometric product of two vectors has a scalar component (determined by their cross product) and a bivector component (determined by their wedge product - an operation similar to the cross product but returns an area spanned by two vectors - this generalizes more easily to higher dimensions than a cross product).

Therefore a quaternion can be thought of as having a unitless scalar component and its i, j, and k components as areas. The conjugation operation multiplies a vector by the quaternion on one side and its inverse on the other, so the units cancel and only the vector's units remain.

giordano commented 1 year ago

While we're at it, speculating about somewhat breaking changes, I was convinced that Quantity should be limited to Real instead of Number, the reasoning being that for example for complex numbers you'd have Complex{Quantity}, instead of the current Quantity{Complex}: you have a complex number of unitful real components, instead of a unitful complex number itself. All units I know are "real". Would there any preferred interpretation for quaternions?

brainandforce commented 1 year ago

This is a really interesting question, and after thinking about it for a bit, I realized that it actually poses a serious problem for someone wanting to incorporate it into a geometric algebra package. (Quaternions are a subalgebra of the 3D vector geometric algebra, so perhaps this is useful to packages that work with quaternions - and complex numbers are a subalgebra of the 2D vector geometric algebra)

I'm currently working on a package that provides multivectors (elements of geometric algebras), which consist of a linear combination of scalars, vectors, bivectors, etc. (up to D-vectors where D is the dimension of the space). That means that the components of the multivector are not all commensurate in units - for instance, with a multivector with length units, the scalar is unitless, the vector is a length, the bivector is an area, etc.

If I wanted to implement a unitful multivector, the way I'd go about it is to give the entire multivector a single unit (basically have a Quantity{<:AbstractMultivector, <:Dimension, <:Units} and work out the specifics of operations on multivectors with new methods as needed to handle the different unit exponents of each component. While it does make sense for AbstractMultivector to subtype Number (multivectors are sometimes called Clifford numbers, after all), limiting the storage type of AbstractQuantity to Real would make this very difficult or cumbersome, if not impossible.

(After attempting to implement Unitful into a different package, I was actually a little surprised that Quantity only wrapped numeric types, and didn't allow for the wrapping of matrices, for instance. I get that it's a subtype of Number and likely benefits from the mathematical operations defined on Number subtypes, but struct definitions can get complicated...)

I understand if this particular issue is too niche to be a significant factor in design considerations, but if there seems to be multiple different possible implementations, I hope this narrows it down a bit.

sostock commented 1 year ago

for complex numbers you'd have Complex{Quantity}

For that, we would need Quantity <: Real, and Iā€™m not sure we should stretch the definition of real numbers that far.