Is there a \ operator for Unitful?

New to Julia
, ,
1

I was running into a problem using the \ (linear system solver) operator in Julia. I guess I wrongly assumed the \ operator was overloaded for matrices/vectors with entries of type <: Number (witch unitful types are). Is there any way to easily overload the \ operator to accept unitful types, or perhaps a different operator/function that works for all types that are subtypes of Number (since I plan on using Sym and Float/Integer types in my code as well)

Environment: Pluto on MacOS
Packages: Unitful, Sympy, Plots, PlutoUI
Code:

let
	a11 = 1u"m"
	a12 = 2u"ft"
	a21 = 5u"mm"
	a22 = 0u"inch"
	A = [a11 a12; a21 a22]
	b = [3u"inch"; 1u"mm"]
	x  = A \ b # This line throws error
end

Error:
MethodError: no method matching (Unitful.Quantity{Rational{Int64}})(::Rational{Int64})

1 Like
2

I don’t think there is an easy way to define a unique, unambiguous solution to that example problem. Which should the units of the elements of x? Moreover, in this case it’s only a matter of scale between units of the same kind, but for arbitrary units in A and b, it might be impossible to solve.

1 Like
3

Maybe something along these lines:

using LinearAlgebra, Unitful
import LinearAlgebra: \

(\)(A::Matrix{Quantity{T,U,V}}, B::AbstractVecOrMat) where {T,U,V} =
  Quantity.(ustrip(A)\ustrip(B), unit(eltype(B))/unit(eltype(A)))

With this definition:

julia> A\b
2-element Vector{Rational{Int64}}:
    1//5
 -619//3048

julia> c = ustrip(b)
2-element reinterpret(Rational{Int64}...):
 381//5000
   1//1000

julia> A\c
2-element Vector{Quantity{Rational{Int64}...}}:
       1//5 m^-1
 -619//3048 m^-1

But this was a quick sketch, might be ins-and-outs to fix yet.

4

Oh yeah, it makes sense that the ‘b’ vector would be in units of Length*mysteryunit since Ax = b, and x is unknown. Maybe x can be solved given the units of b?

5

Using the one-argument ustrip method is generally a bad idea: you don’t know what you’re stripping, it works well only for the trivial case where you have all homogeneous quantities.

2 Likes
6

Thanks for the advice! I changed up the method a little bit to ensure that all unit of elements in the matrices A and B are the same before using the system solver.

function Base.:\(A::Matrix{Quantity{T,U,V}}, B::AbstractVecOrMat) where {T,U,V}
		if typeof(unit(eltype(A))) ≠ Unitful.FreeUnits{(), NoDims, nothing}
			A = map(x -> uconvert(unit(eltype(A)), x), A) 
			B = map(x -> uconvert(unit(eltype(B)), x), B)
			return Quantity.(ustrip(A)\ustrip(B), unit(eltype(B))/unit(eltype(A)))
		end
	end
7

There should be another definition of \ covering the case when A is a non-Unitful matrix and B is Unitful. Should be pretty similar to the one already defined.