ggebbie / UnitfulLinearAlgebra.jl

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

Bug report: `UnitfulMatrix(x)` where `x` is not a `Vector` causes `(Matrix)(::Nothing)` error #92

Closed singularitti closed 1 year ago

singularitti commented 1 year ago

Code example:

julia> using UnitfulLinearAlgebra, StaticArrays

julia> A = SMatrix{3,3}([2 0 0
           0 2 0
           0 0 2] * u"m"
       )
3Γ—3 SMatrix{3, 3, Quantity{Int64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}, 9} with indices SOneTo(3)Γ—SOneTo(3):
 2 m  0 m  0 m
 0 m  2 m  0 m
 0 m  0 m  2 m

julia> x = A * [1, 2, 3]
3-element SVector{3, Quantity{Int64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}} with indices SOneTo(3):
 2 m
 4 m
 6 m

julia> UnitfulMatrix(A) \ UnitfulMatrix(x)
ERROR: MethodError: no method matching (Matrix)(::Nothing)

Closest candidates are:
  (Matrix)(::UndefInitializer, ::Integer, ::Integer)
   @ Base baseext.jl:33
  (Matrix)(::LinearAlgebra.LDLt)
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.0-beta3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/ldlt.jl:223
  (Matrix)(::SparseArrays.AbstractSparseMatrixCSC{Tv}) where Tv
   @ SparseArrays ~/.julia/juliaup/julia-1.10.0-beta3+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/SparseArrays/src/sparsematrix.jl:956
  ...

Stacktrace:
 [1] UnitfulMatrix(a::SVector{3, Quantity{Int64, 𝐋, Unitful.FreeUnits{(m,), 𝐋, nothing}}})
   @ UnitfulLinearAlgebra ~/.julia/packages/UnitfulLinearAlgebra/Rj3QW/src/UnitfulMatrix.jl:172

The cause is because of this line:

if eltype(dims) <: Vector
    return UnitfulMatrix(data, format(Units.(dims), data), exact)
elseif eltype(dims) <: Units
    return UnitfulMatrix(data, format(dims, data), exact)
end

The reason is that you check against eltype(dims) is a subtype of Vector, not AbstractVector, which caused the checking procedure to return a nothing since nothing is matched.

You may want to add an else to notify unexpected behavior:

if eltype(dims) <: Vector
    return UnitfulMatrix(data, format(Units.(dims), data), exact)
elseif eltype(dims) <: Units
    return UnitfulMatrix(data, format(dims, data), exact)
else
    error("this should never happen!")
end

Why is eltype(dims) not a Vector? While, it was generated by this code:

unitrange = Vector{Unitful.Units}(undef,M)

unitrange = unit.(a)

As you can see, unitrange was generated twice, while the first statement is overwritten by the second, for a, i.e., x::SVector{3}, this will cause unitrange to be a SVector, too. Which is obviously not a subtype of Vector.

julia> unit.(a)
3-element SVector{3, Unitful.FreeUnits{(m,), 𝐋, nothing}} with indices SOneTo(3):
 m
 m
 m

You may want to fix it by:

unitrange = Vector{Unitful.Units}(undef,M)

unitrange[:] = unit.(a)
ggebbie commented 1 year ago

Thank you for the detailed synopsis of the situation. I appreciate it!

Note that there may be a whole class of bugs related to StaticArrays that I have not tested. I hope not, but I haven't done any tests except on regular matrices and SparseArrays.

Note that it would be nice to include tests with StaticArrays in runtests but I don't know how to do so without adding the StaticArrays dependency. (In the reverse, I am testing matrices made with SparseArrays, but that dependency could probably be dropped.)

PR #93 is great, but requires three small changes to get the runtests to pass. The biggest issue is that I needed FreeUnits in L169 of UnitfulMatrix.jl. This is a temporary solution, because there are other types of units that are not handled. Issue #88 is probably related.

I tried and failed to update PR #93 (I have not done this before) and instead I have incorporated #93 into #94 which passes tests and passes the StaticArrays code above.