Merge pull request #107 from ggebbie/update-dd

Dimensional Data v0.29

Project.toml

@@ -1,11 +1,12 @@
 name = "UnitfulLinearAlgebra"
 uuid = "c14bd059-d406-4571-8f61-9bd20e53c30b"
 authors = ["G Jake Gebbie <[email protected]>"]
-version = "0.3.8"
+version = "0.4.0"
 
 [deps]
 DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
@@ -17,7 +18,9 @@ UnitfulLatexify = "45397f5d-5981-4c77-b2b3-fc36d6e9b728"
 UnitfulLinearAlgebraLatexifyExt = ["Latexify", "UnitfulLatexify"]
 
 [compat]
-DimensionalData = "0.24,0.25,0.27"
+#DimensionalData = "0.24,0.25,0.27,0.28"
+DimensionalData = "0.29"
+Revise = "3.6.2"
 Statistics = "1"
 Unitful = "1"
 julia = "1.9,1.10,1.11"

src/UnitfulLinearAlgebra.jl

@@ -27,7 +27,7 @@ export describe, trace
 import Base: (\)
 import LinearAlgebra: inv
 import DimensionalData: @dim, dims, DimArray, AbstractDimArray, NoName, NoMetadata, format, print_name
-#import Latexify: latexify
+#import Unitful: Units #import Latexify: latexify
 
 @dim Units "units"
 

src/UnitfulMatrix.jl

@@ -92,13 +92,17 @@ update is dealt with in `rebuild` for `AbstractDimArray` (still true?).
     urange = unitrange(A)[I[1]]
     udomain = unitdomain(A)[I[2]]
 
-    if (udomain isa Unitful.Units || urange isa Unitful.Units )
-        # case of column vector, row vector, scalar
-        # scalar appears to be overridden by getindex
-        newunitrange = slicedvector(urange,udomain)
+    # something strange with two types of `Units` type
+    if (udomain isa Unitful.Units && urange isa UnitfulLinearAlgebra.Units )
+        # case of column vector
+        newunitrange = slicedvector(parent(urange),udomain)
+        return UnitfulMatrix(data, newunitrange,exact=false)
+    elseif (udomain isa UnitfulLinearAlgebra.Units && urange isa Unitful.Units )
+        # case of row vector
+        newunitrange = slicedvector(urange,parent(udomain))
         return UnitfulMatrix(data, newunitrange,exact=false)
-        #return UnitfulMatrix(data, newunitrange, newunitdomain)
     else
+        # matrix case
         newunitrange, newunitdomain = slicedmatrix(urange,udomain)
         # unit range and domain of a sliced matrix are ambiguous.
         # It must be exact=false

src/base.jl

@@ -19,7 +19,6 @@ function Base.show(io::IO, mime::MIME"text/plain", B::AbstractUnitfulDimVecOrMat
     return nothing
 end
 
-
 function Base.show(io::IO, mime::MIME"text/plain", A::AbstractUnitfulVecOrMat)
     lines = 0
     summary(io, A)
@@ -75,9 +74,23 @@ Base.:*(b::Number,a::Union{AbstractUnitfulVector,AbstractUnitfulDimVector}) = a*
 # (matrix/vector)-(matrix/vector) multiplication when inexact handled here
 function Base.:*(A::AbstractUnitfulVecOrMat,B::AbstractUnitfulVecOrMat)
     if exact(A) && exact(B)
-        return DimensionalData._rebuildmul(A,B)
+        # replace DimensionalData._rebuildmul(A,B) # uses strict checking
+        # instead reproduce necessary part of DimensionalData here
+
+        # from DimensionalData._comparedims_mul(A, B)
+        DimensionalData.comparedims(last(dims(A)), first(dims(B)); 
+            order=false, val=true, length=false
+        )
+        return rebuild(A, parent(A) * parent(B), (first(dims(A)), last(dims(B))))
+
     elseif unitdomain(A) ∥ unitrange(B)
-        return DimensionalData._rebuildmul(convert_unitdomain(A,unitrange(B)),B)
+
+        Anew = convert_unitdomain(A,unitrange(B))
+        DimensionalData.comparedims(last(dims(Anew)),
+            first(dims(B)); 
+            order=false, val=true, length=false
+        )
+        return rebuild(Anew, parent(Anew) * parent(B), (first(dims(Anew)), last(dims(B))))
     else
         error("UnitfulLinearAlgebra: unitdomain(A) and unitrange(B) not parallel")
     end
@@ -124,7 +137,7 @@ Base.:-(A::AbstractUnitfulType) = DimensionalData.rebuild(A,-parent(A))
      Reverse mapping from unitdomain to range.
      Is `exact` if input is exact.
 """
-function Base.:\(A::AbstractUnitfulMatrix,b::AbstractUnitfulVector)
+function Base.:(\ )(A::AbstractUnitfulMatrix,b::AbstractUnitfulVector)
     if exact(A)
         DimensionalData.comparedims(first(dims(A)), first(dims(b)); val=true)
         return rebuild(A,parent(A)\parent(b),(last(dims(A)),)) #,exact = (exact(A) && exact(B)))
@@ -135,7 +148,7 @@ function Base.:\(A::AbstractUnitfulMatrix,b::AbstractUnitfulVector)
         error("UnitfulLinearAlgebra.mldivide: Dimensions of Unitful Matrices A and b not compatible")
     end
 end
-function Base.:\(A::AbstractUnitfulMatrix,B::AbstractUnitfulMatrix)
+function Base.:(\ )(A::AbstractUnitfulMatrix,B::AbstractUnitfulMatrix)
     if exact(A)
         DimensionalData.comparedims(first(dims(A)), first(dims(B)); val=true)
         return rebuild(A,parent(A)\parent(B),(last(dims(A)),last(dims(B)))) #,exact = (exact(A) && exact(B)))
@@ -146,7 +159,7 @@ function Base.:\(A::AbstractUnitfulMatrix,B::AbstractUnitfulMatrix)
         error("UnitfulLinearAlgebra.matrix left divide): Dimensions of Unitful Matrices A and b not compatible")
     end
 end
-function Base.:\(A::AbstractUnitfulDimMatrix,b::AbstractUnitfulDimVector)
+function Base.:(\ )(A::AbstractUnitfulDimMatrix,b::AbstractUnitfulDimVector)
     if exact(A)
         DimensionalData.comparedims(first(unitdims(A)), first(unitdims(b)); val=true)
         DimensionalData.comparedims(first(dims(A)), first(dims(b)); val=true)
@@ -158,7 +171,7 @@ function Base.:\(A::AbstractUnitfulDimMatrix,b::AbstractUnitfulDimVector)
         error("UnitfulLinearAlgebra.mldivide: Dimensions of Unitful Matrices A and b not compatible")
     end
 end
-function Base.:\(A::AbstractUnitfulDimMatrix,B::AbstractUnitfulDimMatrix)
+function Base.:(\ )(A::AbstractUnitfulDimMatrix,B::AbstractUnitfulDimMatrix)
     if exact(A)
         DimensionalData.comparedims(first(unitdims(A)), first(unitdims(B)); val=true)
         DimensionalData.comparedims(first(dims(A)), first(dims(B)); val=true)
@@ -171,10 +184,8 @@ function Base.:\(A::AbstractUnitfulDimMatrix,B::AbstractUnitfulDimMatrix)
     end
 end
 # do what the investigator means -- convert to UnitfulType -- probably a promotion mechanism to do the same thing
-Base.:\(A::AbstractUnitfulType,b::Number) = A\UnitfulMatrix([b])
-# this next one is quite an assumption
-#Base.:\(A::AbstractUnitfulMatrix,b::AbstractVector) = A\UnitfulMatrix(vec(b)) #error("UnitfulLinearAlgebra: types not consistent")
-Base.:\(A::AbstractUnitfulMatrix,b::Vector) = vec(A\UnitfulMatrix(b)) # return something with same type as input `b`
+Base.:(\ )(A::AbstractUnitfulType,b::Number) = A\UnitfulMatrix([b])
+Base.:(\ )(A::AbstractUnitfulMatrix,b::Vector) = vec(A\UnitfulMatrix(b)) # return something with same type as input `b`
 
 """
     function ldiv(F::LU{T,MultipliableMatrix{T},Vector{Int64}}, B::AbstractVector) where T<:Number
@@ -207,41 +218,6 @@ function (\)(F::LU{T,<: AbstractUnitfulMatrix,Vector{Int64}}, B::AbstractUnitful
     return rebuild(B,LinearAlgebra._cut_B(BB, 1:n),(unitdomain(F.factors),))
 end
 
-# """
-#      function ldiv!
-
-#      In-place left division by a Multipliable Matrix.
-#      Reverse mapping from unitdomain to range.
-#      Is `exact` if input is exact.
-
-#     Problem: b changes type unless endomorphic
-# """
-# function ldiv!(A::AbstractMultipliableMatrix,b::AbstractVector)
-#     ~endomorphic(A) && error("A not endomorphic, b changes type, ldiv! not available")
-
-#     if dimension(unitrange(A)) == dimension(b)
-#         #if unitrange(A) ~ b
-
-#         # seems to go against the point
-#         #b = copy((A.numbers\ustrip.(b)).*unitdomain(A))
-#         btmp = (A.numbers\ustrip.(b)).*unitdomain(A)
-#         for bb = 1:length(btmp)
-#             b[bb] = btmp[bb]
-#         end
-
-#     elseif ~exact(A) && (unitrange(A) ∥ b)
-#         Anew = convert_unitrange(A,unit.(b)) # inefficient?
-#         btmp = (Anew.numbers\ustrip.(b)).*unitdomain(Anew)
-#         for bb = 1:length(btmp)
-#             b[bb] = btmp[bb]
-#         end
-
-#     else
-#         error("UnitfulLinearAlgebra.ldiv!: Dimensions of MultipliableMatrix and vector not compatible")
-#     end
-
-# end
-
 Base.:~(a,b) = similarity(a,b)
 
 """
@@ -252,8 +228,9 @@ Base.:~(a,b) = similarity(a,b)
 
     Hart, pp. 205.
 """
-# A redefined tranpose that corrects error based on AbstractArray interface
-Base.transpose(A::AbstractUnitfulMatrix) = rebuild(A,transpose(parent(A)),(Units(unitdomain(A).^-1), Units(unitrange(A).^-1)))
+# previous working version here
+#Base.transpose(A::AbstractUnitfulMatrix) = rebuild(A,transpose(parent(A)),(Units(inv.(unitdomain(A))), Units(inv.(unitrange(A)))))
+Base.transpose(A::AbstractUnitfulMatrix) = rebuild(A,transpose(parent(A)),(Units(parent(inv.(unitdomain(A)))), Units(parent(inv.(unitrange(A))))))
 Base.transpose(a::AbstractUnitfulVector) = rebuild(a,transpose(parent(a)),(Units([NoUnits]), Units(unitrange(a).^-1))) # kludge for unitrange of row vector
 Base.transpose(A::AbstractUnitfulDimMatrix) = rebuild(A,transpose(parent(A)),(Units(unitdomain(A).^-1), Units(unitrange(A).^-1)),(last(dims(A)),first(dims(A))))
 Base.transpose(a::AbstractUnitfulDimVector) = rebuild(a,transpose(parent(a)),(Units([NoUnits]), Units(unitrange(a).^-1)),(:empty,first(dims(a))))
@@ -271,34 +248,6 @@ Base.similar(A::AbstractUnitfulVecOrMat{T}) where T <: Number =
 
 # NOTE: Base.getproperty is also expanded but stored next to relevant linear algebra functions.
 
-# """
-#     function vcat(A,B)
-
-#     Modeled after function `VERTICAL` (pp. 203, Hart, 1995).
-# """
-# function Base.vcat(A::AbstractMultipliableMatrix,B::AbstractMultipliableMatrix)
-
-#     numbers = vcat(A.numbers,B.numbers)
-#     shift = unitdomain(A)[1]./unitdomain(B)[1]
-#     ur = vcat(unitrange(A),unitrange(B).*shift)
-#     bothexact = (exact(A) && exact(B))
-#     return BestMultipliableMatrix(numbers,ur,unitdomain(A),exact=bothexact)
-# end
-
-# """
-#     function hcat(A,B)
-
-#     Modeled after function `HORIZONTAL` (pp. 202, Hart, 1995).
-# """
-# function Base.hcat(A::AbstractMultipliableMatrix,B::AbstractMultipliableMatrix)
-
-#     numbers = hcat(A.numbers,B.numbers)
-#     shift = unitrange(A)[1]./unitrange(B)[1]
-#     ud = vcat(unitdomain(A),unitdomain(B).*shift)
-#     bothexact = (exact(A) && exact(B))
-#     return BestMultipliableMatrix(numbers,unitrange(A),ud,exact=bothexact)
-# end
-
 ## start of UnitfulDimMatrix methods
 Base.:*(A::AbstractUnitfulDimMatrix, B::AbstractUnitfulDimMatrix) = DimensionalData._rebuildmul(A,B)
 
@@ -396,3 +345,32 @@ end
 
 Base.first(A::AbstractUnitfulType) = first(vec(A))
 Base.last(A::AbstractUnitfulType) = last(vec(A))
+
+
+# """
+#     function vcat(A,B)
+
+#     Modeled after function `VERTICAL` (pp. 203, Hart, 1995).
+# """
+# function Base.vcat(A::AbstractMultipliableMatrix,B::AbstractMultipliableMatrix)
+
+#     numbers = vcat(A.numbers,B.numbers)
+#     shift = unitdomain(A)[1]./unitdomain(B)[1]
+#     ur = vcat(unitrange(A),unitrange(B).*shift)
+#     bothexact = (exact(A) && exact(B))
+#     return BestMultipliableMatrix(numbers,ur,unitdomain(A),exact=bothexact)
+# end
+
+# """
+#     function hcat(A,B)
+
+#     Modeled after function `HORIZONTAL` (pp. 202, Hart, 1995).
+# """
+# function Base.hcat(A::AbstractMultipliableMatrix,B::AbstractMultipliableMatrix)
+
+#     numbers = hcat(A.numbers,B.numbers)
+#     shift = unitrange(A)[1]./unitrange(B)[1]
+#     ud = vcat(unitdomain(A),unitdomain(B).*shift)
+#     bothexact = (exact(A) && exact(B))
+#     return BestMultipliableMatrix(numbers,unitrange(A),ud,exact=bothexact)
+# end

src/dsvd.jl

@@ -83,7 +83,7 @@ Base.propertynames(F::DSVD, private::Bool=false) =
 
     Dimensioned singular value decomposition (DSVD).
     Appropriate version of SVD for non-uniform matrices.
-    `svd` can be computed for `Number`s, `Adjoint`s, `Tranpose`s, and `Integers`; `dsvd` doesn't yet implement these.
+    `svd` can be computed for `Number`s, `Adjoint`s, `Transpose`s, and `Integers`; `dsvd` doesn't yet implement these.
 # Input
 - `A::AbstractMultipliableMatrix`
 - `Pr::UnitSymmetricMatrix`: square matrix defining norm of range
@@ -107,13 +107,8 @@ function dsvd(A::AbstractUnitfulMatrix,Py::AbstractUnitfulMatrix,Px::AbstractUni
 
     # matrix slices cause unit domain and range to become ambiguous.
     # output of DSVD cannot be exact.
-    # if exact(A)
-    #     U = convert_unitdomain(UnitfulMatrix(F.U),Units(fill(unit(1.0),size(F.U,2))))
-    #     Vt = convert_unitrange(UnitfulMatrix(F.Vt),Units(fill(unit(1.0),size(F.Vt,1))))
-    #     return DSVD(U,F.S,Vt,Qy,Qx)
-    # else
-        return DSVD(UnitfulMatrix(F.U),F.S,UnitfulMatrix(F.Vt),Qy,Qx)
-    #end
+    return DSVD(UnitfulMatrix(F.U),F.S,UnitfulMatrix(F.Vt),Qy,Qx)
+
 end
 
 function show(io::IO, mime::MIME{Symbol("text/plain")}, F::DSVD{<:Any,<:Any,<:AbstractArray,<:AbstractArray,<:AbstractArray,<:AbstractArray,<:AbstractVector})
@@ -133,30 +128,18 @@ function inv(F::DSVD{T}) where T
         iszero(F.S[i]) && throw(SingularException(i))
     end
     k = searchsortedlast(F.S, eps(real(T))*F.S[1], rev=true)
-    # adapted from `svd.jl`
-    #@views (F.S[1:k] .\ F.V[:,1:k, :]) #* F.U⁻¹[1:k,:]
 
+    # adapted from `svd.jl`
     # make sure it is inexact
     Σ⁻¹ = UnitfulMatrix(Diagonal(F.S[1:k].^-1))
-    # a less efficient matrix way to do it.
-    #Σ⁻¹ = Diagonal(F.S[1:k].^-1,fill(unit(1.0),k),fill(unit(1.0),k))
-    # Σ⁻¹ = Diagonal(F.S[1:k].^-1,unitdomain(F.V[:,1:k]),unitrange(F.U⁻¹[1:k,:]))
-    #    Σ⁻¹ = Diagonal(F.S[1:k].^-1,unitdomain(F.V)[1:k],unitrange(F.U⁻¹)[1:k])
-    #println(exact(F.V[:,1:k]))
     return F.V[:,1:k]*Σ⁻¹*F.U⁻¹[1:k,:]
 end
 
-### DSVD least squares ### Not implemented
-# function ldiv!(A::SVD{T}, B::StridedVecOrMat) where T
-#     m, n = size(A)
-#     k = searchsortedlast(A.S, eps(real(T))*A.S[1], rev=true)
-#     mul!(view(B, 1:n, :), view(A.Vt, 1:k, :)', view(A.S, 1:k) .\ (view(A.U, :, 1:k)' * _cut_B(B, 1:m)))
-#     return B
-# end
-
 Base.size(A::DSVD, dim::Integer) = dim == 1 ? size(A.U, dim) : size(A.V⁻¹, dim)
 Base.size(A::DSVD) = (size(A, 1), size(A, 2))
 
+### DSVD least squares ### Not implemented
+
 # adjoint not yet defined for AbstractMultipliableMatrix
 #function adjoint(F::DSVD)
 #    return SVD(F.V⁻¹', F.S, F.U')

src/linear_algebra.jl

@@ -10,8 +10,6 @@
 """
 LinearAlgebra.inv(A::AbstractUnitfulMatrix) = rebuild(A,inv(parent(A)),(unitdomain(A),unitrange(A)))
 LinearAlgebra.inv(A::AbstractUnitfulDimMatrix) = rebuild(A,inv(parent(A)), (unitdomain(A),unitrange(A)), (last(dims(A)),first(dims(A)) ))
-# PR 78: singular check sometimes fails but matrix is invertible.
-#LinearAlgebra.inv(A::AbstractUnitfulMatrix) = ~singular(A) ? rebuild(A,inv(parent(A)),(unitdomain(A),unitrange(A))) : error("matrix is singular")
 
 """
     function det
@@ -68,10 +66,10 @@ function LinearAlgebra.inv(A::Eigen{T,V,S,U}) where {U <: AbstractVector, S <: A
         ur = unitrange(A.vectors)
         ud = Units(unit.(A.values))
         Λ⁻¹ = Diagonal(A.values.^-1,ud,ur)
-        return A.vectors* transpose(transpose(A.vectors) \ Λ⁻¹)
+        #return A.vectors* transpose(transpose(A.vectors) \ Λ⁻¹)
 
         # LinearAlgebra.eigen uses matrix right divide, i.e., 
-        #return A.vectors * Λ⁻¹ / A.vectors
+        return A.vectors * Λ⁻¹ / A.vectors
         # but this is not available yet for `Multipliable Matrix`s.
 
     else
@@ -107,7 +105,8 @@ end
 function LinearAlgebra.cholesky(A::AbstractUnitfulMatrix)
     if unit_symmetric(A)
         C = LinearAlgebra.cholesky(parent(A))
-        factors = rebuild(A,C.factors,(Units(unitdomain(A)./unitdomain(A)),unitdomain(A)))
+        #factors = rebuild(A,C.factors,(Units(unitdomain(A)./unitdomain(A)),unitdomain(A)))
+        factors = rebuild(A,C.factors,(Units(fill(NoUnits,size(A,1))),unitdomain(A)))
         return Cholesky(factors,C.uplo,C.info)
     else
         error("requires unit symmetric matrix")
@@ -116,17 +115,14 @@ end
 function LinearAlgebra.cholesky(A::AbstractUnitfulDimMatrix)
     if unit_symmetric(A)
         C = LinearAlgebra.cholesky(parent(A))
-
-        # What should the axis units for a Cholesky decomposition be? Just a guess here.
-        # What if internal parts of Cholesky decomposition are simply UnitfulMatrix's. 
         factors = rebuild(A,C.factors,(Units(unitdomain(A)./unitdomain(A)),unitdomain(A)),(:Normalspace,last(dims(A))))
         return Cholesky(factors,C.uplo,C.info)
     else
         error("requires unit symmetric matrix")
     end
 end
 
-# seems like this might be from Base? Move to base.jl? 
+# Move to base.jl? 
 function Base.getproperty(C::Cholesky{T,<:AbstractUnitfulMatrix}, d::Symbol) where T 
     Cfactors = getfield(C, :factors)
     Cuplo    = getfield(C, :uplo)