Unitful and generic function types

General Usage
, ,
1

I have a function not using Unitful that reads like

function foo(x :: Matrix{T}, y :: Vector{T}) where T
    my_array = zeros(T, 100)
    # other stuff 
end

Now, I want to implement my code using Unitful to make it easier to spot bugs with dimensional checks, allow the user to pass data in whatever units he/she likes the most, and have Unitful handle the conversion.
However, I am struggling to get information about the underlying type. For example, it is nice that the function checks that the type of elements of x and y are the same when called. And I also can use the T to initialize new vectors inside the function to the correct type.

However, not necessary the x, y might have the same units, in that case, Unitful will return an error. I can avoid the check, but then also initializing my_array to something with different units from x and y is a nightmare.
I have to do something like (if I want an array of seconds)

my_array = zeros(typeof(zero(eltype(ustrip(x))) * u"s"), 100)

Which is definitively horrible to read and understand.

Is there a proper way to manage the numeric type besides the units and easily create and manage arrays with Unitful?

2

One option is to use AbstractMatrix{T} and AbstractVector{R} rather than using T for both. That way you can get the type without requiring the arguments have the same units.

Another option would be to convert both arguments to the same units before calling this function.

3

It’s a little weird to me that ustrip doesn’t take the Quantity type as an input like unit can, so this would be what that would look like:

julia> extra_ustrip(::Type{Q}) where {T, Q<:Quantity{T}} = T
extra_ustrip (generic function with 1 method)

julia> extra_ustrip(typeof(1*u"s")) # replace input with `T` in your code
Int64

julia> extra_ustrip(typeof(1im*u"s"))
Complex{Int64}

But I don’t know how that would be multiplied with u"s".

You don’t provide much detail, so let’s say we write the results of bar(x, y) to my_array. For that to work, we need eltype(my_array) == eltype(bar(x,y)). This isn’t easy to deal with automatically; the output type of the first elementwise operation is not necessarily the same as the following ones, which is why the output type mechanisms for broadcasting are complex and stayed internal. The common practice is to pass the output eltype or the preallocated output array itself to foo, leaving another function to figure out the automatic output type mechanism for your own context or allowing you to manually provide the expected output type for x and y.

4

Thanks @Benny and @jar1
Looking at your code, probably the solution is replacing my definition of foo into

function foo(x :: AbstractMatrix{T}, y :: AbstractVector{U}) where {V, T <: Quantity{V}, U <: Quantity{V}}
    unit_type = typeof(zero(V) * u"s")
    my_array = zeros(unit_type, 100)
    # other stuff
end

In this way, I am forcing the underlying float representation of input to be the same. I wonder if this function signature is generic enough to work with GPU arrays.