Tidy one() implementation and UniformScaling constructors (#690)

* Remove all reference to the deprecated name `eye` (replace internal
  _eye with _scaler_matrix)
* Remove all the more-or-less duplicate definitions of one() and use
  _scalar_matrix as the implementation.
* Fix StaticMatrixLike definition to properly include Hermitian and
  Symmetric wrappers of static matrices.
* Import LinearAlgebra.checksquare() and use it directly.
* Remove obsolete checksquare specializations & add tests
* Introduce _construct_similar to alleviate the difficulty of computing
  eltype purely in the type domain with non-concrete UnionAlls.

Author: c42f

src/MArray.jl

@@ -73,16 +73,6 @@ end
 
 @inline MArray(a::StaticArray) = MArray{size_tuple(Size(a))}(Tuple(a))
 
-# Simplified show for the type
-#show(io::IO, ::Type{MArray{S, T, N}}) where {S, T, N} = print(io, "MArray{$S,$T,$N}")
-
-# Some more advanced constructor-like functions
-@inline one(::Type{MArray{S}}) where {S} = one(MArray{S,Float64,tuple_length(S)})
-@inline one(::Type{MArray{S,T}}) where {S,T} = one(MArray{S,T,tuple_length(S)})
-
-# MArray(I::UniformScaling) methods to replace eye
-(::Type{MA})(I::UniformScaling) where {MA<:MArray} = _eye(Size(MA), MA, I)
-
 ####################
 ## MArray methods ##
 ####################

src/MVector.jl

@@ -95,7 +95,7 @@ macro MVector(ex)
                 error("@MVector expected a 1-dimensional array expression")
             end
         else
-            error("@MVector only supports the zeros(), ones(), rand(), randn(), randexp(), and eye() functions.")
+            error("@MVector only supports the zeros(), ones(), rand(), randn(), and randexp() functions.")
         end
     else
         error("Use @MVector [a,b,c] or @MVector([a,b,c])")

src/SArray.jl

@@ -52,16 +52,6 @@ end
 
 @inline SArray(a::StaticArray) = SArray{size_tuple(Size(a))}(Tuple(a))
 
-# Simplified show for the type
-# show(io::IO, ::Type{SArray{S, T, N}}) where {S, T, N} = print(io, "SArray{$S,$T,$N}") # TODO reinstate
-
-# Some more advanced constructor-like functions
-@inline one(::Type{SArray{S}}) where {S} = one(SArray{S, Float64, tuple_length(S)})
-@inline one(::Type{SArray{S, T}}) where {S, T} = one(SArray{S, T, tuple_length(S)})
-
-# SArray(I::UniformScaling) methods to replace eye
-(::Type{SA})(I::UniformScaling) where {SA<:SArray} = _eye(Size(SA), SA, I)
-
 ####################
 ## SArray methods ##
 ####################

src/SDiagonal.jl

@@ -33,7 +33,7 @@ size(::Type{<:SDiagonal{N}}, d::Int) where {N} = d > 2 ? 1 : N
 # override to avoid copying
 diag(D::SDiagonal) = D.diag
 
-# SDiagonal(I::UniformScaling) methods to replace eye
+# SDiagonal(I::UniformScaling) methods
 (::Type{SDiagonal{N}})(I::UniformScaling) where {N} = SDiagonal{N}(ntuple(x->I.λ, Val(N)))
 (::Type{SDiagonal{N,T}})(I::UniformScaling) where {N,T} = SDiagonal{N,T}(ntuple(x->I.λ, Val(N)))
 

src/SHermitianCompact.jl

@@ -230,7 +230,7 @@ end
     end
 end
 
-@inline _eye(s::Size{S}, t::Type{SSC}) where {S, SSC <: SHermitianCompact} = _one(s, t)
+@inline _scalar_matrix(s::Size{S}, t::Type{SSC}) where {S, SSC <: SHermitianCompact} = _one(s, t)
 
 # _fill covers fill, zeros, and ones:
 @generated function _fill(val, ::Size{s}, ::Type{SSC}) where {s, SSC <: SHermitianCompact}

src/SVector.jl

@@ -108,9 +108,9 @@ macro SVector(ex)
                 error("@SVector expected a 1-dimensional array expression")
             end
         else
-            error("@SVector only supports the zeros(), ones(), rand(), randn(), randexp(), and eye() functions.")
+            error("@SVector only supports the zeros(), ones(), rand(), randn() and randexp() functions.")
         end
-    else # TODO Expr(:call, :zeros), Expr(:call, :ones), Expr(:call, :eye) ?
+    else
         error("Use @SVector [a,b,c], @SVector Type[a,b,c] or a comprehension like [f(i) for i = i_min:i_max]")
     end
 end

src/StaticArrays.jl

@@ -19,11 +19,7 @@ import LinearAlgebra: transpose, adjoint, dot, eigvals, eigen, lyap, tr,
                       kron, diag, norm, dot, diagm, lu, svd, svdvals,
                       factorize, ishermitian, issymmetric, isposdef, normalize,
                       normalize!, Eigen, det, logdet, cross, diff, qr, \
-
-# import eye for deprecation warnings
-@static if isdefined(LinearAlgebra, :eye)
-    import LinearAlgebra: eye
-end
+using LinearAlgebra: checksquare
 
 export SOneTo
 export StaticScalar, StaticArray, StaticVector, StaticMatrix
@@ -88,8 +84,8 @@ const StaticMatrixLike{T} = Union{
     StaticMatrix{<:Any, <:Any, T},
     Transpose{T, <:StaticVecOrMat{T}},
     Adjoint{T, <:StaticVecOrMat{T}},
-    Symmetric{T, <:StaticMatrix{T}},
-    Hermitian{T, <:StaticMatrix{T}},
+    Symmetric{T, <:StaticMatrix{<:Any, <:Any, T}},
+    Hermitian{T, <:StaticMatrix{<:Any, <:Any, T}},
     Diagonal{T, <:StaticVector{<:Any, T}}
 }
 const StaticVecOrMatLike{T} = Union{StaticVector{<:Any, T}, StaticMatrixLike{T}}

src/abstractarray.jl

@@ -93,6 +93,27 @@ mutable_similar_type(::Type{T}, s::Size{S}, ::Type{Val{D}}) where {T,S,D} = MArr
 
 sizedarray_similar_type(::Type{T},s::Size{S},::Type{Val{D}}) where {T,S,D} = SizedArray{Tuple{S...},T,D,length(s)}
 
+# Utility for computing the eltype of an array instance, type, or type
+# constructor.  For type constructors without a definite eltype, the default
+# value is returned.
+Base.@pure _eltype_or(a::AbstractArray, default) = eltype(a)
+Base.@pure _eltype_or(::Type{<:AbstractArray{T}}, default) where {T} = T
+Base.@pure _eltype_or(::Type{<:AbstractArray}, default) = default # eltype not available
+
+"""
+    _construct_similar(a, ::Size, elements::NTuple)
+
+Construct a static array of similar type to `a` with the given `elements`.
+
+When `a` is an instance or a concrete type the element type `eltype(a)` is
+used. However, when `a` is a `UnionAll` type such as `SMatrix{2,2}`, the
+promoted type of `elements` is used instead.
+"""
+@inline function _construct_similar(a, s::Size, elements::NTuple{L,ET}) where {L,ET}
+    similar_type(a, _eltype_or(a, ET), s)(elements)
+end
+
+
 # Field vectors are user controlled, and currently default to SVector, etc
 
 """

src/det.jl

@@ -43,7 +43,7 @@ end
 end
 
 @generated function _det(::Size{S}, A::StaticMatrix) where S
-    LinearAlgebra.checksquare(A)
+    checksquare(A)
     if prod(S) ≤ 14*14
         quote
             @_inline_meta
@@ -58,7 +58,7 @@ end
 @inline logdet(A::StaticMatrix) = _logdet(Size(A), A)
 @inline _logdet(::Union{Size{(1,1)}, Size{(2,2)}, Size{(3,3)}, Size{(4,4)}}, A::StaticMatrix) = log(det(A))
 @generated function _logdet(::Size{S}, A::StaticMatrix) where S
-    LinearAlgebra.checksquare(A)
+    checksquare(A)
     if prod(S) ≤ 14*14
         quote
             @_inline_meta

src/linalg.jl

@@ -1,6 +1,7 @@
 import Base: +, -, *, /, \
 
-# TODO: more operators, like AbstractArray
+#--------------------------------------------------
+# Vector space algebra
 
 # Unary ops
 @inline -(a::StaticArray) = map(-, a)
@@ -30,10 +31,7 @@ import Base: +, -, *, /, \
 @inline -(a::UniformScaling, b::StaticMatrix) = _plus_uniform(Size(b), -b, a.λ)
 
 @generated function _plus_uniform(::Size{S}, a::StaticMatrix, λ) where {S}
-    if S[1] != S[2]
-        throw(DimensionMismatch("matrix is not square: dimensions are $S"))
-    end
-    n = S[1]
+    n = checksquare(a)
     exprs = [i == j ? :(a[$(LinearIndices(S)[i, j])] + λ) : :(a[$(LinearIndices(S)[i, j])]) for i = 1:n, j = 1:n]
     return quote
         $(Expr(:meta, :inline))
@@ -46,6 +44,8 @@ end
 @inline \(a::UniformScaling, b::Union{StaticMatrix,StaticVector}) = a.λ \ b
 @inline /(a::StaticMatrix, b::UniformScaling) = a / b.λ
 
+#--------------------------------------------------
+# Matrix algebra
 
 # Transpose, conjugate, etc
 @inline conj(a::StaticArray) = map(conj, a)
@@ -85,31 +85,30 @@ end
 @inline Base.zero(a::SA) where {SA <: StaticArray} = zeros(SA)
 @inline Base.zero(a::Type{SA}) where {SA <: StaticArray} = zeros(SA)
 
-@inline one(::SM) where {SM <: StaticMatrix} = _one(Size(SM), SM)
-@inline one(::Type{SM}) where {SM <: StaticMatrix} = _one(Size(SM), SM)
-@generated function _one(::Size{S}, ::Type{SM}) where {S, SM <: StaticArray}
-    if (length(S) != 2) || (S[1] != S[2])
-        error("multiplicative identity defined only for square matrices")
-    end
-    T = eltype(SM) # should be "hyperpure"
-    if T == Any
-        T = Float64
-    end
-    exprs = [i == j ? :(one($T)) : :(zero($T)) for i ∈ 1:S[1], j ∈ 1:S[2]]
-    return quote
-        $(Expr(:meta, :inline))
-        SM(tuple($(exprs...)))
+@inline one(m::StaticMatrixLike) = _one(Size(m), m)
+@inline one(::Type{SM}) where {SM<:StaticMatrixLike}= _one(Size(SM), SM)
+function _one(s::Size, m_or_SM)
+    if (length(s) != 2) || (s[1] != s[2])
+        throw(DimensionMismatch("multiplicative identity defined only for square matrices"))
     end
+    _scalar_matrix(s, m_or_SM, one(_eltype_or(m_or_SM, Float64)))
 end
 
-# StaticMatrix(I::UniformScaling) methods to replace eye
-(::Type{SM})(I::UniformScaling) where {N,M,SM<:StaticMatrix{N,M}} = _eye(Size(SM), SM, I)
-
-@generated function _eye(::Size{S}, ::Type{SM}, I::UniformScaling{T}) where {S, SM <: StaticArray, T}
-    exprs = [i == j ? :(I.λ) : :(zero($T)) for i ∈ 1:S[1], j ∈ 1:S[2]]
+# StaticMatrix(I::UniformScaling)
+(::Type{SM})(I::UniformScaling) where {SM<:StaticMatrix} = _scalar_matrix(Size(SM), SM, I.λ)
+# The following oddity is needed if we want `SArray{Tuple{2,3}}(I)` to work
+# because we do not have `SArray{Tuple{2,3}} <: StaticMatrix`.
+(::Type{SM})(I::UniformScaling) where {SM<:(StaticArray{Tuple{N,M}} where {N,M})} =
+    _scalar_matrix(Size(SM), SM, I.λ)
+
+# Construct a matrix with the scalar λ on the diagonal and zeros off the
+# diagonal. The matrix can be non-square.
+@generated function _scalar_matrix(s::Size{S}, m_or_SM, λ) where {S}
+    elements = Symbol[i == j ? :λ : :λzero for i in 1:S[1], j in 1:S[2]]
     return quote
         $(Expr(:meta, :inline))
-        SM(tuple($(exprs...)))
+        λzero = zero(λ)
+        _construct_similar(m_or_SM, s, tuple($(elements...)))
     end
 end
 
@@ -145,6 +144,8 @@ end
     end
 end
 
+#--------------------------------------------------
+# Vector products
 @inline cross(a::StaticVector, b::StaticVector) = _cross(same_size(a, b), a, b)
 _cross(::Size{S}, a::StaticVector, b::StaticVector) where {S} = error("Cross product not defined for $(S[1])-vectors")
 @inline function _cross(::Size{(2,)}, a::StaticVector, b::StaticVector)
@@ -179,6 +180,8 @@ end
     return ret
 end
 
+#--------------------------------------------------
+# Norms
 @inline LinearAlgebra.norm_sqr(v::StaticVector) = mapreduce(abs2, +, v; init=zero(real(eltype(v))))
 
 @inline norm(a::StaticArray) = _norm(Size(a), a)
@@ -240,9 +243,7 @@ end
 
 @inline tr(a::StaticMatrix) = _tr(Size(a), a)
 @generated function _tr(::Size{S}, a::StaticMatrix) where {S}
-    if S[1] != S[2]
-        throw(DimensionMismatch("matrix is not square"))
-    end
+    checksquare(a)
 
     if S[1] == 0
         return :(zero(eltype(a)))
@@ -257,6 +258,10 @@ end
     end
 end
 
+
+#--------------------------------------------------
+# Outer products
+
 const _length_limit = Length(200)
 
 @inline kron(a::StaticMatrix, b::StaticMatrix) = _kron(_length_limit, Size(a), Size(b), a, b)
@@ -414,11 +419,9 @@ end
     end
 end
 
-# some micro-optimizations (TODO check these make sense for v0.6+)
-@inline LinearAlgebra.checksquare(::SM) where {SM<:StaticMatrix} = _checksquare(Size(SM))
-@inline LinearAlgebra.checksquare(::Type{SM}) where {SM<:StaticMatrix} = _checksquare(Size(SM))
 
-@pure _checksquare(::Size{S}) where {S} = (S[1] == S[2] || throw(DimensionMismatch("matrix is not square: dimensions are $S")); S[1])
+#--------------------------------------------------
+# Some shimming for special linear algebra matrix types
+@inline LinearAlgebra.Symmetric(A::StaticMatrix, uplo::Char='U') = (checksquare(A); Symmetric{eltype(A),typeof(A)}(A, uplo))
+@inline LinearAlgebra.Hermitian(A::StaticMatrix, uplo::Char='U') = (checksquare(A); Hermitian{eltype(A),typeof(A)}(A, uplo))
 
-@inline LinearAlgebra.Symmetric(A::StaticMatrix, uplo::Char='U') = (LinearAlgebra.checksquare(A);Symmetric{eltype(A),typeof(A)}(A, uplo))
-@inline LinearAlgebra.Hermitian(A::StaticMatrix, uplo::Char='U') = (LinearAlgebra.checksquare(A);Hermitian{eltype(A),typeof(A)}(A, uplo))