Inconsistency with ustrip(array) and mixed units

[sostock]

When calling ustrip on a quantity with mixed units like dBV/m, both logarithmic and linear units are stripped. When calling ustrip on an array of these quantities, only the linear part is stripped:

julia> x = 5u"dBV/m"
[5.0 dBV] m^-1

julia> ustrip(x)
5.0

julia> ustrip([x])
1-element reinterpret(Level{Unitful.LogInfo{:Decibel, 10, 10}, 1 V, Quantity{Float64, 𝐋^2 𝐌 𝐈^-1 𝐓^-3, Unitful.FreeUnits{(V,), 𝐋^2 𝐌 𝐈^-1 𝐓^-3, nothing}}}, ::Vector{Quantity{Level{Unitful.LogInfo{:Decibel, 10, 10}, 1 V, Quantity{Float64, 𝐋^2 𝐌 𝐈^-1 𝐓^-3, Unitful.FreeUnits{(V,), 𝐋^2 𝐌 𝐈^-1 𝐓^-3, nothing}}}, 𝐋 𝐌 𝐈^-1 𝐓^-3, Unitful.FreeUnits{(m^-1,), 𝐋^-1, nothing}}}):
 5.0 dBV

This is because calling ustrip on an Array{<:Quantity} actually reinterprets the array. It is not possible to reinterpret the value 5.0dBV/m to the value 5.0, because the Level type stores the linear quantity (in this case 1.7782794100389228 V, not the number 5.0).

[bjarthur]

ustrip is nominally deprecated, for reasons i don't understand. try your example using broadcasting instead:

julia> x = 5u"dBV/m"
[5.0 dBV] m^-1

julia> ustrip.(x)
5.0

julia> ustrip.([x])
1-element Vector{Float64}:
 5.0

the downside of broadcasting is that ustrip then allocates:

julia> aa=rand(1000)*1u"ms";

julia> @time ustrip(aa);
  0.000010 seconds (1 allocation: 32 bytes)

julia> @time ustrip.(aa);
  0.000022 seconds (3 allocations: 7.969 KiB)

[bjarthur]

wait, hold up! the method of ustrip that is deprecated is the one that inputs a non-unitful array. still, i wonder why, because broadcasting ustrip onto a non-unitful array allocates, and it'd be nice to avoid that when writing generic code.

[sostock]

It’s not great that ustrip(::Array) reinterprets in some cases and allocates a new array in other cases. The option to reinterpret is useful, for example if you want to pass the array to a BLAS function.

I propose the following for Unitful 2.0:

• If you don’t want to reinterpret, just use broadcasting (this already works).
• If you want to reinterpret, you should specify it explicitly (to avoid surprises). Something like ustrip(reinterpret, x) looks good to me. So we would remove all ustrip(::SomeArray) methods and add ustrip(::typeof(reinterpret), ::SomeArray) methods instead. These would error if the correct value cannot be obtained by reinterpreting (e.g. for mixed units).

[bjarthur]

by SomeArray do you mean to include AbstractArray? i ask, because it sure would be nice if Unitful played better with CuArrays. for example, currently, one cannot ustrip a CuArray and then assign to it:

julia> using Unitful, CUDA

julia> c = CuArray(zeros(3) * 1u"s")
3-element CuArray{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}, 1, CUDA.Mem.DeviceBuffer}:
 0.0 s
 0.0 s
 0.0 s

julia> ustrip(c)[1]=1   # innocuous warning about scalar indexing removed for succinctness
1

julia> c
3-element CuArray{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}, 1, CUDA.Mem.DeviceBuffer}:
 0.0 s
 0.0 s
 0.0 s

however, if i simply add a new method identical to this one but with CuArray, then it works fine:

julia> import Unitful: ustrip

julia> @inline ustrip(A::CuArray{Q}) where {Q <: Quantity} = reinterpret(Unitful.numtype(Q), A)
ustrip (generic function with 16 methods)

julia> ustrip(c)[1]=1
1

julia> c
3-element CuArray{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}, 1, CUDA.Mem.DeviceBuffer}:
 1.0 s
 0.0 s
 0.0 s

to be all inclusive, i think the right way to go here is just use AbstractArray. then it'll work too for MtlArray, ROCArray, etc.

[bjarthur]

here's another pain point with composing Unitful and CUDA. randn gives different outputs if the input array has units but is on a CPU vs a GPU.

everything is fine w/o Unitful. randn's output is the same given the same seed:

julia> using CUDA, Random

julia> len=30;

julia> rng = MersenneTwister();

julia> a = zeros(len);

julia> Random.seed!(rng, 1);

julia> randn!(rng, a)
30-element Vector{Float64}:
  0.2972879845354616
  0.3823959677906078
 -0.5976344767282311
 -0.01044524463737564
 -0.839026854388764
  0.31111133849833383
  2.2950878238373105
 -0.05317733752769253
  0.2290095549097807
 -2.2670863488005306
  0.5299655761667461
  0.43142152642291204
  0.5837082875687786
  0.9632716050381906
  0.45879095505371686
 -0.5223367574215084
  0.40839583832475224
 -0.050451229933665284
 -0.6936536438038856
 -1.7738332534080556
  0.12076596493743928
 -0.7576330377205277
 -1.7297561813906863
  0.7959486220046159
  0.6700619812560624
  0.5508523775007609
 -0.06337459242956062
  1.3369427822509223
 -0.07314863333773934
 -0.7454643166407556

julia> rng
MersenneTwister(1, (0, 1002, 0, 32))

julia> c = CuArray(zeros(len));

julia> Random.seed!(rng, 1);

julia> randn!(rng, c)
30-element CuArray{Float64, 1, CUDA.Mem.DeviceBuffer}:
  0.2972879845354616
  0.3823959677906078
 -0.5976344767282311
 -0.01044524463737564
 -0.839026854388764
  0.31111133849833383
  2.2950878238373105
 -0.05317733752769253
  0.2290095549097807
 -2.2670863488005306
  0.5299655761667461
  0.43142152642291204
  0.5837082875687786
  0.9632716050381906
  0.45879095505371686
 -0.5223367574215084
  0.40839583832475224
 -0.050451229933665284
 -0.6936536438038856
 -1.7738332534080556
  0.12076596493743928
 -0.7576330377205277
 -1.7297561813906863
  0.7959486220046159
  0.6700619812560624
  0.5508523775007609
 -0.06337459242956062
  1.3369427822509223
 -0.07314863333773934
 -0.7454643166407556

julia> a==c   # innocuous warning removed again there
true

julia> rng
MersenneTwister(1, (0, 1002, 0, 32))

doing the same thing as above, but this time with both arrays having units, and using the ustrip method in my post above to reinterpret so that an assignment is possible, results in different random numbers. curiously, the first seven are the same. and the state of rng afterwards is also the same.

julia> using Unitful

julia> import Unitful: ustrip

julia> @inline ustrip(A::CuArray{Q}) where {Q <: Quantity} = reinterpret(Unitful.numtype(Q), A)
ustrip (generic function with 16 methods)

julia> ua = zeros(len) * 1u"s";

julia> Random.seed!(rng, 1);

julia> randn!(rng, ustrip(ua))
30-element reinterpret(Float64, ::Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}):
  0.2972879845354616
  0.3823959677906078
 -0.5976344767282311
 -0.01044524463737564
 -0.839026854388764
  0.31111133849833383
  2.2950878238373105
 -2.2670863488005306
  0.5299655761667461
  0.43142152642291204
  0.5837082875687786
  0.9632716050381906
  0.45879095505371686
 -0.5223367574215084
  0.40839583832475224
 -0.050451229933665284
 -0.6936536438038856
 -1.7738332534080556
  0.12076596493743928
 -0.7576330377205277
 -1.7297561813906863
  0.7959486220046159
  0.6700619812560624
  0.5508523775007609
 -0.06337459242956062
  1.3369427822509223
 -0.07314863333773934
 -0.7454643166407556
 -1.2200551285346526
 -0.05317733752769253

julia> ua
30-element Vector{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}}:
    0.2972879845354616 s
    0.3823959677906078 s
   -0.5976344767282311 s
  -0.01044524463737564 s
    -0.839026854388764 s
   0.31111133849833383 s
    2.2950878238373105 s
   -2.2670863488005306 s
    0.5299655761667461 s
   0.43142152642291204 s
    0.5837082875687786 s
    0.9632716050381906 s
   0.45879095505371686 s
   -0.5223367574215084 s
   0.40839583832475224 s
 -0.050451229933665284 s
   -0.6936536438038856 s
   -1.7738332534080556 s
   0.12076596493743928 s
   -0.7576330377205277 s
   -1.7297561813906863 s
    0.7959486220046159 s
    0.6700619812560624 s
    0.5508523775007609 s
  -0.06337459242956062 s
    1.3369427822509223 s
  -0.07314863333773934 s
   -0.7454643166407556 s
   -1.2200551285346526 s
  -0.05317733752769253 s

julia> rng
MersenneTwister(1, (0, 1002, 0, 32))

julia> uc = CuArray(zeros(len) * 1u"s");

julia> Random.seed!(rng, 1);

julia> randn!(rng, ustrip(uc))
30-element CuArray{Float64, 1, CUDA.Mem.DeviceBuffer}:
  0.2972879845354616
  0.3823959677906078
 -0.5976344767282311
 -0.01044524463737564
 -0.839026854388764
  0.31111133849833383
  2.2950878238373105
 -0.05317733752769253
  0.2290095549097807
 -2.2670863488005306
  0.5299655761667461
  0.43142152642291204
  0.5837082875687786
  0.9632716050381906
  0.45879095505371686
 -0.5223367574215084
  0.40839583832475224
 -0.050451229933665284
 -0.6936536438038856
 -1.7738332534080556
  0.12076596493743928
 -0.7576330377205277
 -1.7297561813906863
  0.7959486220046159
  0.6700619812560624
  0.5508523775007609
 -0.06337459242956062
  1.3369427822509223
 -0.07314863333773934
 -0.7454643166407556

julia> uc
30-element CuArray{Quantity{Float64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}, 1, CUDA.Mem.DeviceBuffer}:
    0.2972879845354616 s
    0.3823959677906078 s
   -0.5976344767282311 s
  -0.01044524463737564 s
    -0.839026854388764 s
   0.31111133849833383 s
    2.2950878238373105 s
  -0.05317733752769253 s
    0.2290095549097807 s
   -2.2670863488005306 s
    0.5299655761667461 s
   0.43142152642291204 s
    0.5837082875687786 s
    0.9632716050381906 s
   0.45879095505371686 s
   -0.5223367574215084 s
   0.40839583832475224 s
 -0.050451229933665284 s
   -0.6936536438038856 s
   -1.7738332534080556 s
   0.12076596493743928 s
   -0.7576330377205277 s
   -1.7297561813906863 s
    0.7959486220046159 s
    0.6700619812560624 s
    0.5508523775007609 s
  -0.06337459242956062 s
    1.3369427822509223 s
  -0.07314863333773934 s
   -0.7454643166407556 s

julia> ua==uc
false

julia> rng
MersenneTwister(1, (0, 1002, 0, 32))

[sostock]

by SomeArray do you mean to include AbstractArray?

I was mainly thinking about the structured matrices for which we already have specialized methods:

Unitful.jl/src/utils.jl

Lines 89 to 92 in 69ffe80

	ustrip(A::Diagonal) = Diagonal(ustrip(A.diag))
	ustrip(A::Bidiagonal) = Bidiagonal(ustrip(A.dv), ustrip(A.ev), ifelse(istriu(A), :U, :L))
	ustrip(A::Tridiagonal) = Tridiagonal(ustrip(A.dl), ustrip(A.d), ustrip(A.du))
	ustrip(A::SymTridiagonal) = SymTridiagonal(ustrip(A.dv), ustrip(A.ev))

But it would be nice to define it for generic AbstractArrays as well.

I don’t know whether defining a generic AbstractArray method is enough to make it work correctly for CuArrays (I haven’t used CUDA and don’t own an NVIDIA GPU, so I cannot test it).

[bjarthur]

It’s not great that ustrip(::Array) reinterprets in some cases and allocates a new array in other cases.

the only case where ustrip(::Array) allocates that i can find is a method that was deprecated in v0.4. we are hence free to delete and/or redefine that definition without violating semvar.

without that method, ustrip on an array always reinterprets so far as i can tell. if the user would rather allocate, they should be instructed to broadcast instead.

i'm not sure i see the need to add ustrip(reinterpret... methods. we should rather just delete that deprecation to make the interface consistent.

[sostock]

The fact that the method is deprecated does not mean that we can remove it without making a breaking release. Deprecating just means “we recommend not using this method and may remove it in the next breaking release”, so we would still have to release v2.0 if we want to remove that method.

The reason why I would like to introduce ustrip(reinterpret, x) is that you can easily mistake ustrip(x) for ustrip.(x). Getting an informative error when you forget the dot seems better to me than getting unexpected results by mistakenly reinterpreting.

[bjarthur]

my bad about semvar deprecations. you're right.

re. missing a dot to broadcast-- i don't think many julian developers would make that mistake. it's too common in the language.