FFT for Unitful vector

General Usage
,
1

I am trying to use Unitful package to better track units of my physical quantities, and I am a bit confused:

Here is example:

using Unitful

trajectory = trajectory .* u"m"

fft(trajectory)

I got the following error:

ERROR: MethodError: no method matching plan_fft(::Matrix{Quantity{Float64, ๐‹, Unitful.FreeUnits{(m,), ๐‹, nothing}}}, ::UnitRange{Int64})

It is a bug or a feature?

2

To use quantities with units, packages need to be written in a generic way. You didnโ€™t say where youโ€™re getting the function fft from, but if itโ€™s from say FFTW, then thatโ€™s calling out to a library written in C, and that code can only handle a few hardware-supported numeric types like Float64. So you need to unwrap and rewrap your quantities to interface with such code, e.g. fft(ustrip.(u"m", trajectory)).

You can actually still get benefits from using units in this way though, e.g.

julia> x = rand(4)
4-element Vector{Float64}:
 0.16106228971764325
 0.7034731999527042
 0.6401812191825158
 0.43054032170468615

julia> x_m = x * u"m"
4-element Vector{Quantity{Float64, ๐‹, Unitful.FreeUnits{(m,), ๐‹, nothing}}}:
 0.16106228971764325 m
  0.7034731999527042 m
  0.6401812191825158 m
 0.43054032170468615 m

julia> x_mm = x * u"mm"
4-element Vector{Quantity{Float64, ๐‹, Unitful.FreeUnits{(mm,), ๐‹, nothing}}}:
 0.16106228971764325 mm
  0.7034731999527042 mm
  0.6401812191825158 mm
 0.43054032170468615 mm

julia> ustrip.(u"m", x_m)
4-element Vector{Float64}:
 0.16106228971764325
 0.7034731999527042
 0.6401812191825158
 0.43054032170468615

julia> ustrip.(u"m", x_mm)
4-element Vector{Float64}:
 0.00016106228971764325
 0.0007034731999527042
 0.0006401812191825159
 0.00043054032170468615

That is, ustrip will convert quantities to the type you specify before stripping off the unit, so that you can be sure what goes into fft has the units you expect.

5 Likes
3

Notice BTW that you can strip (and replace) units without copying โ€“ they are purely in the type:

julia> x_re = reinterpret(Float64, x_m)
4-element reinterpret(Float64, ::Vector{Quantity{Float64, ๐‹, Unitful.FreeUnits{(m,), ๐‹, nothing}}}):
 0.16106228971764325
 0.7034731999527042
 0.6401812191825158
 0.43054032170468615

julia> x_re[1] # no units, but a view of x_m
0.16106228971764325

julia> fft(x_re) # accepts ReinterpretArray
4-element Vector{ComplexF64}:
   1.9352570305575494 + 0.0im
 -0.47911892946487256 - 0.27293287824801804im
  -0.3327700127572313 + 0.0im
 -0.47911892946487256 + 0.27293287824801804im

julia> @which fft
AbstractFFTs

Also, the basic definitions of fft etc. are in a small package to make it easy to extend them, without the package which does so being tied to FFTW.

7 Likes
4

Probably works only for arrays of homogeneous quantities (all with exactly the same units)?

5 Likes
5

Yes, absolutely. I guess it could still make sense to take the FFT of something like Any[1u"m", 1u"cm"], for which this trick would not work.

(And, I think fft copies real vectors to make complex ones before getting to work, in any case. So ideally one would do both steps at once, and reinterpret only dense arrays of unit-homogeneous complex numbers.)

1 Like
6

Thank you all for answers!

How should I deal with matrix where are NaNs?

a = zeros(N,1) .*u"m*s^-2"
if somecondition
  a[ii] = value_in_ms^-2
else
  a[ii] = NaN
end
7

You would probably have to impute it. For instance, have its value entirely depend on its neighbors (e.g. with Interpolations.jl).

8

Thatโ€™s super useful! I was considering whether or not use to use units in a package where I wasnโ€™t sure if the cost of allocating a new array was worth it, but this makes it basically free.

Would it make sense to put something like this in Unitful itself?

"""
    with_unit(v::AbstractArray, u::Unitful.Units) -> Base.ReinterpretArray{โ€ฆ}

Creates a view of `v` in which each element has units given by `u`.

!!! warning
    This does not perform any unit conversion or check the units of `v` first.
"""
with_unit(v::AbstractArray, u::Unitful.Units) = reinterpret(typeof(one(eltype(v))*u), v)
9

You could even spell it x_m = reinterpret(u"m", x) so as not to need to remember a new name; perhaps a slight abuse but not piracy.

Should any existing functions default to this (when the array is suitable)? Unitful doesnโ€™t seem to have much support for arrays, things like uconvert.(u"cm", x_m) and ustrip.(u"cm", x_m) are at present errors without the broadcasting.

10

There are a few linear algebra operations that donโ€™t work with unitful arrays:

And stripping the units isnโ€™t always a good idea, especially for inhomogenous arrays (another example in the issue)

11

I am struggling with the type of argument in the function

I would like to construct multiple dispatch function

for example:

>> a = rand(N) .*u"m*s^-2"

>> typeof(a)
Vector{Quantity{Float64, ๐‹ ๐“^-2, Unitful.FreeUnits{(m, s^-2), ๐‹ ๐“^-2, nothing}}}

function do_st_with_a(a)
 do_st_with_a # units are not needed
end

function do_st_with_a(a:: Vector{Quantity{Float64, ๐‹ ๐“^-2, Unitful.FreeUnits{(m, s^-2), ๐‹ ๐“^-2, nothing}}})
 do_st_with_a # units are needed
end

This give me an error: ERROR: syntax: missing comma or } in argument list

What type should I put to the second definition of the function to get multiple dispatch depending on presence of units in input array?

The second thing which confused me is this

>> typeof(a[1]) == typeof(Unitful.Acceleration)

false