Merge branch 'JuliaArrays:master' into master

Author: astrozot

Project.toml

@@ -1,6 +1,6 @@
 name = "StaticArrays"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.8"
+version = "1.9.13"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

ext/StaticArraysStatisticsExt.jl

@@ -1,9 +1,12 @@
 module StaticArraysStatisticsExt
 
-import Statistics: mean
+import Statistics: mean, median
+
+using Base.Order: Forward, ord
+using Statistics: median!, middle
 
 using StaticArrays
-using StaticArrays: _InitialValue, _reduce, _mapreduce
+using StaticArrays: BitonicSort, _InitialValue, _reduce, _mapreduce, _bitonic_sort_limit, _sort
 
 _mean_denom(a, ::Colon) = length(a)
 _mean_denom(a, dims::Int) = size(a, dims)
@@ -12,4 +15,37 @@ _mean_denom(a, ::Val{D}) where {D} = size(a, D)
 @inline mean(a::StaticArray; dims=:) = _reduce(+, a, dims) / _mean_denom(a, dims)
 @inline mean(f::Function, a::StaticArray; dims=:) = _mapreduce(f, +, dims, _InitialValue(), Size(a), a) / _mean_denom(a, dims)
 
+@inline function median(a::StaticArray; dims = :)
+    if dims == Colon()
+        median(vec(a))
+    else
+        # FIXME: Implement `mapslices` correctly on `StaticArray` to remove
+        # this fallback.
+        median(Array(a); dims)
+    end
+end
+
+@inline function median(a::StaticVector)
+    (isimmutable(a) && length(a) <= _bitonic_sort_limit) ||
+        return median!(Base.copymutable(a))
+
+    # following Statistics.median
+    isempty(a) &&
+        throw(ArgumentError("median of empty vector is undefined, $(repr(a))"))
+    eltype(a) >: Missing && any(ismissing, a) &&
+        return missing
+    nanix = findfirst(x -> x isa Number && isnan(x), a)
+    isnothing(nanix) ||
+        return a[nanix]
+
+    order = ord(isless, identity, nothing, Forward)
+    sa = _sort(Tuple(a), BitonicSort, order)
+
+    n = length(a)
+    # sa is 1-indexed
+    return isodd(n) ?
+        middle(sa[n ÷ 2 + 1]) :
+        middle(sa[n ÷ 2], sa[n ÷ 2 + 1])
+end
+
 end # module

src/MArray.jl

@@ -42,7 +42,7 @@ end
     return v
 end
 
-@inline Tuple(v::MArray) = getfield(v,:data)
+@inline Base.Tuple(v::MArray) = getfield(v,:data)
 
 Base.dataids(ma::MArray) = (UInt(pointer(ma)),)
 

src/MVector.jl

@@ -18,6 +18,11 @@ end
 let dimension_names = QuoteNode.([:x, :y, :z, :w])
     body = :(getfield(v, name))
     for (i,dim_name) in enumerate(dimension_names)
+        @eval @inline function Base.propertynames(v::Union{SVector{$i},MVector{$i}}, private::Bool = false)
+            named_dims = ($(first(dimension_names, i)...),)
+            private ? (named_dims..., :data) : named_dims
+        end
+
         body = :(name === $(dimension_names[i]) ? getfield(v, :data)[$i] : $body)
         @eval @inline function Base.getproperty(v::Union{SVector{$i},MVector{$i}},
                                                 name::Symbol)
@@ -33,3 +38,6 @@ let dimension_names = QuoteNode.([:x, :y, :z, :w])
         end
     end
 end
+
+# for longer S/MVectors and other S/MArrays, the only property is data, and it's private
+Base.propertynames(::Union{SArray,MArray}, private::Bool = false) = private ? (:data,) : ()

src/SArray.jl

@@ -65,7 +65,7 @@ sacollect
     getfield(v,:data)[i]
 end
 
-@inline Tuple(v::SArray) = getfield(v,:data)
+@inline Base.Tuple(v::SArray) = getfield(v,:data)
 
 Base.dataids(::SArray) = ()
 

src/convert.jl

@@ -198,7 +198,7 @@ AbstractArray{T}(sa::StaticArray{S,U}) where {S,T,U} = similar_type(typeof(sa),T
 AbstractArray{T,N}(sa::StaticArray{S,U,N}) where {S,T,U,N} = similar_type(typeof(sa),T,Size(sa))(sa)
 
 # Constructing a Tuple from a StaticArray
-@inline Tuple(a::StaticArray) = unroll_tuple(a, Length(a))
+@inline Base.Tuple(a::StaticArray) = unroll_tuple(a, Length(a))
 
 @noinline function dimension_mismatch_fail(::Type{SA}, a::AbstractArray) where {SA <: StaticArray}
     throw(DimensionMismatch("expected input array of length $(length(SA)), got length $(length(a))"))

src/eigen.jl

@@ -400,7 +400,8 @@ end
 @inline eigen(A::Symmetric{<:Complex,<:StaticMatrix}; kwargs...) = _eigen(A; kwargs...)
 
 @inline function _eigen(A::LinearAlgebra.HermOrSym; permute::Bool=true, scale::Bool=true)
-    _eig(Size(A), A, permute, scale)
+    B = convert(AbstractArray{float(eltype(A))}, A)
+    _eig(Size(A), B, permute, scale)
 end
 
 # NOTE: The following Boost Software License applies to parts of the method:

src/expm.jl

@@ -14,6 +14,7 @@ end
     (newtype)((exp(A[1]), ))
 end
 
+# Bernstein, D. S. & So, W. 1993. "Some Explicit Formulas for the Matrix Exponential"
 @inline function _exp(::Size{(2,2)}, A::StaticMatrix{<:Any,<:Any,<:Real})
     T = typeof(exp(zero(eltype(A))))
     newtype = similar_type(A,T)
@@ -25,25 +26,50 @@ end
 
     v = (a-d)^2 + 4*b*c
 
+    m = (a + d) / 2
+
     if v > 0
-        z = sqrt(v)
-        z1 = cosh(z / 2)
-        z2 = sinh(z / 2) / z
+        # In this case the formulas in the entries of the matrix
+        # are a function of cosh and sinh, and could (in theory)
+        # follow the same code pattern of the other branches (v ≤ 0).
+        # However, cosh and sinh explode with large arguments and
+        # we use the following identity to avoid numerical issues:
+        #
+        # exp(m) * [c₁ * cohs(δ) + c₂ * sinh(δ)] =
+        #   c₁ * (e₊ + e₋) / 2 + c₂ * (e₊ - e₋) / 2
+        #
+        # where e₊ = exp(m + δ) and e₋ = exp(m - δ).
+        #
+        # See https://github.com/JuliaArrays/StaticArrays.jl/issues/1295
+        δ = sqrt(v) / 2
+        e₊ = exp(m + δ)
+        e₋ = exp(m - δ)
+        e₁ = (e₊ + e₋) / 2
+        e₂ = (e₊ - e₋) / 2
+        c₂ = (a - d) / 2δ
+        m11 = (e₁ + c₂ * e₂)
+        m12 = (b / δ) * e₂
+        m21 = (c / δ) * e₂
+        m22 = (e₁ - c₂ * e₂)
     elseif v < 0
         z = sqrt(-v)
+        r = exp(m)
         z1 = cos(z / 2)
         z2 = sin(z / 2) / z
+        m11 = r * (z1 + (a - d) * z2)
+        m12 = r * 2b * z2
+        m21 = r * 2c * z2
+        m22 = r * (z1 - (a - d) * z2)
     else # if v == 0
+        r = exp(m)
         z1 = T(1.0)
         z2 = T(0.5)
+        m11 = r * (z1 + (a - d) * z2)
+        m12 = r * 2b * z2
+        m21 = r * 2c * z2
+        m22 = r * (z1 - (a - d) * z2)
     end
 
-    r = exp((a + d) / 2)
-    m11 = r * (z1 + (a - d) * z2)
-    m12 = r * 2b * z2
-    m21 = r * 2c * z2
-    m22 = r * (z1 - (a - d) * z2)
-
     (newtype)((m11, m21, m12, m22))
 end
 

src/linalg.jl

@@ -301,7 +301,7 @@ end
         $(Expr(:meta, :inline))
         scale = maxabs_nested(a)
 
-        scale==0 && return _init_zero(a)
+        iszero(scale) && return _init_zero(a)
         p == 1 && return @inbounds scale * $expr_p1
         return @inbounds scale * ($expr)^(inv(p))
     end
@@ -328,7 +328,7 @@ end
         p == Inf && return mapreduce(norm, max, a)  # no need for scaling
 
         l = p==1 ? @inbounds($expr_p1) : @inbounds(($expr)^(inv(p)))
-        0<l<Inf && return l
+        zero(l) < l && isfinite(l) && return l
         return _norm_scaled(Size(a), a, p)  # p != 0, 2, Inf
     end
 end

src/matrix_multiply.jl

@@ -114,7 +114,7 @@ end
         b::Union{Transpose{Tb, <:StaticVector}, Adjoint{Tb, <:StaticVector}}) where {sa, sb, Ta, Tb}
     newsize = (sa[1], sb[2])
     exprs = [:(a[$i]*b[$j]) for i = 1:sa[1], j = 1:sb[2]]
-    
+
     return quote
         @_inline_meta
         T = promote_op(*, Ta, Tb)
@@ -208,7 +208,7 @@ end
         while m < M
             mu = min(M, m + M_r)
             mrange = m+1:mu
-            
+
             atemps_init = [:($(atemps[k1]) = a[$k1]) for k1 = mrange]
             exprs_init = [:($(tmps[k1,k2])  = $(atemps[k1]) * b[$(1 + (k2-1) * sb[1])]) for k1 = mrange, k2 = nrange]
             atemps_loop_init = [:($(atemps[k1]) = a[$(k1-sa[1]) + $(sa[1])*j]) for k1 = mrange]
@@ -246,7 +246,7 @@ end
         T = promote_op(matprod, Ta, Tb)
         a = mul_parent(wrapped_a)
         b = mul_parent(wrapped_b)
-        return (mul_result_structure(wrapped_a, wrapped_b))(similar_type(a, T, $S)(invoke(*, Tuple{$(_unstatic_array(a)),$(_unstatic_array(b))}, a, b)))
+        return (mul_result_structure(wrapped_a, wrapped_b))(similar_type(a, T, $S)(invoke(*, Tuple{$_unstatic_array(a),$_unstatic_array(b)}, a, b)))
     end
 end