I ran into something interesting yesterday - I wrote this function to create the cross-product matrix for 3-vector (sometimes known as the tilde operator):
tilde(x) = [0 -x[3] x[2]; x[3] 0 -x[1]; -x[2] x[1] 0]
If I then use it with an input that has units, it seemingly works OK:
julia> v = [0.1 0.2 0.3]u"rad/s"
1×3 Array{Quantity{Float64,𝐓⁻¹,Unitful.FreeUnits{(rad, s⁻¹),𝐓⁻¹,nothing}},2}:
0.1 rad s⁻¹ 0.2 rad s⁻¹ 0.3 rad s⁻¹
julia> ṽ = tilde(v)
3×3 Array{Quantity{Float64,D,U} where U where D,2}:
0.0 -0.3 rad s⁻¹ 0.2 rad s⁻¹
0.3 rad s⁻¹ 0.0 -0.1 rad s⁻¹
-0.2 rad s⁻¹ 0.1 rad s⁻¹ 0.0
Trouble begins, though, when I try to use this matrix:
julia> ṽ*ṽ'
ERROR: DimensionError: 0.0 and 0.09 rad² s⁻² are not dimensionally compatible.
Stacktrace:
[1] +(::Quantity{Float64,NoDims,Unitful.FreeUnits{(),NoDims,nothing}}, ::Quantity{Float64,𝐓⁻²,Unitful.FreeUnits{(rad², s⁻²),𝐓⁻²,nothing}}) at /Users/patrick/.julia/packages/Unitful/1t88N/src/quantities.jl:137
[2] +(::Float64, ::Quantity{Float64,𝐓⁻²,Unitful.FreeUnits{(rad², s⁻²),𝐓⁻²,nothing}}) at ./promotion.jl:311
[3] +(::Float64, ::Quantity{Float64,𝐓⁻²,Unitful.FreeUnits{(rad², s⁻²),𝐓⁻²,nothing}}, ::Quantity{Float64,𝐓⁻²,Unitful.FreeUnits{(rad², s⁻²),𝐓⁻²,nothing}}) at ./operators.jl:538
[4] matmul3x3!(::Array{Any,2}, ::Char, ::Char, ::Array{Quantity{Float64,D,U} where U where D,2}, ::Array{Quantity{Float64,D,U} where U where D,2}, ::LinearAlgebra.MulAddMul{true,true,Bool,Bool}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/matmul.jl:996
[5] generic_matmatmul!(::Array{Any,2}, ::Char, ::Char, ::Array{Quantity{Float64,D,U} where U where D,2}, ::Array{Quantity{Float64,D,U} where U where D,2}, ::LinearAlgebra.MulAddMul{true,true,Bool,Bool}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/matmul.jl:726
[6] mul! at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/matmul.jl:411 [inlined]
[7] mul! at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/matmul.jl:208 [inlined]
[8] *(::Array{Quantity{Float64,D,U} where U where D,2}, ::LinearAlgebra.Adjoint{Any,Array{Quantity{Float64,D,U} where U where D,2}}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/matmul.jl:153
[9] top-level scope at REPL[10]:100:
I didn’t want to write the function in a way that requires Unitful.jl and found that this does the trick:
function tilde(x)
z = zero(x[1]) # trick to make this play nice with Unitful.jl
[z -x[3] x[2]; x[3] z -x[1]; -x[2] x[1] z]
end
The output of this has the same units for every entry and so multiplies without error.
I guess the general question this raises to me is whether there’s a best practice in general for writing a function like mine so that it will play nice with Unitful (and possibly other packages that extend quantities in Julia in interesting ways).
This also got me thinking about whether there’s a way for a vector or matrix as a whole to have units rather than each entry individually having its own units. This does indeed seem to be the case as the type of the array is different when every entry has the same units (or if the units are specified for the array as a whole):
julia> [1u"m" 1u"s"] # entries with different units
1×2 Array{Quantity{Int64,D,U} where U where D,2}:
1 m 1 s
julia> [1u"m" 1u"m"] # all entries same units
1×2 Array{Quantity{Int64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},2}:
1 m 1 m
julia> [1 1]u"m" # entries unit-less, array given units
1×2 Array{Quantity{Int64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},2}:
1 m 1 m
I’m not sure this can help in my example given the way the matrix is constructed – In the function, how could Julia know that what I want is a matrix with all entries the same type as those of x? Only, I guess, if the code that constructs such a matrix were to assume any zero entry has the same units as the nonzero entries, whenever the latter have the same units for each entry.
Incidentally I wondered if the two ways of creating the “homogeneous units” array were the same under the hood; surprisingly (or maybe not) specifying units on each entry rather than for the whole array seems significantly faster:
julia> fa() = [1u"m" 2u"m"]
fa (generic function with 1 method)
julia> fb() = [1 2]u"m"
fb (generic function with 1 method)
julia> fa() == fb()
true
julia> @btime fa()
29.271 ns (1 allocation: 96 bytes)
1×2 Array{Quantity{Int64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},2}:
1 m 2 m
julia> @btime fb()
69.644 ns (2 allocations: 192 bytes)
1×2 Array{Quantity{Int64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}},2}:
1 m 2 m
Also relatedly: Obviously for 3x3 cross product matrices there’s not much performance gain to be had but I’m wondering why there isn’t a skew symmetric special matrix type in the LinearAlgebra package.