ggebbie / UnitfulLinearAlgebra.jl

Low-cost linear algebra functions for matrices with units
MIT License
6 stars 3 forks source link

dsvd doesn't work for FixedUnits #88

Closed b-r-hamilton closed 1 year ago

b-r-hamilton commented 1 year ago

The following code block

using Unitful, UnitfulLinearAlgebra, LinearAlgebra
permil = Unitful.FixedUnits(u"permille")
urange = fill(permil, 2)
udomain = fill(permil, 3)
A = UnitfulMatrix(randn(length(urange),length(udomain)),urange,udomain)
Pr = UnitfulMatrix(I(length(urange)),urange.^-1,urange)
Pd = UnitfulMatrix(I(length(udomain)),udomain.^-1,udomain)
U,Λ,V = dsvd(A,Pr,Pd)

throws the following error

ERROR: DimensionMismatch: Dimension index [, ] and Unitful.Units[, ] do not match

This seems to trace back to the fact that when the following is called https://github.com/ggebbie/UnitfulLinearAlgebra.jl/blob/de6592342d2aaf7c2eaae4d59125cc2b04a8f690/src/dsvd.jl#L61 F.Qy and F.U′ have different unitranges/domains, so the ldiv can't be called. E.g.

julia> svd = dsvd(A,Pr,Pd)
julia> svd.Qy.dims
Units Unitful.FixedUnits{(), NoDims, nothing}[, ],
Units Sampled{Unitful.FixedUnits{(‰,), NoDims, nothing}} Unitful.FixedUnits{(‰,), NoDims, nothing}[‰, ‰] Unordered Irregular Points

julia> svd.U′.dims
Units Categorical{Unitful.Units} Unitful.Units[, ] Unordered,
Units Categorical{Unitful.Units} Unitful.Units[, ] Unordered

It seems to me like the Categorical Units versus the Sampled Units is the problem, and I'm not sure how to fix this. Would there be a way to expand on the test of equivalent dimensions that ldiv calls to make it broader? Like if the vector version of the units is equivalent, it's okay to ldiv?

ggebbie commented 1 year ago

converting the unitrange of Uprime is one solution. dsvd.jl is updated in PR #100 to do this for all cases, which may not be necessary and may slow things down.

There is a similar issue with Vprime in DSVD solved with same approach in PR #100 .

A wider solution could be: 1) keep Uprime and Vprime not UnitfulMatrixs (just Matrix{Quantity} instead), 2) make sure ldiv works properly between a UnitfulMatrix and Matrix{Quantity} (including updating the exact criteria for multiplication, including automatically converting the unitranges of inexact matrices).

b-r-hamilton commented 1 year ago

The code snippet I commented here now runs with the new UnitfulLinearAlgebra. It took a little work to get BLUEs.jl/runtests.jl to work with it, but that's also running with new ULA.