Merge branch 'master' into matrix-exponential-master

Godisemo · web-flow · commit b31a48e3520e · 2017-05-23T09:24:51.000+02:00
diff --git a/src/FieldVector.jl b/src/FieldVector.jl
@@ -23,8 +23,6 @@ abstract type FieldVector{N, T} <: StaticVector{N, T} end
 @propagate_inbounds setindex!(v::FieldVector, x, i::Int) = setfield!(v, i, x)
 
 # See #53
-Base.cconvert{T}(::Type{Ptr{T}}, v::FieldVector) = Base.RefValue(v)
-Base.unsafe_convert{T, FV <: FieldVector}(::Type{Ptr{T}}, m::Base.RefValue{FV}) =
-    _unsafe_convert(Ptr{T}, eltype(FV), m)
-_unsafe_convert{T, FV <: FieldVector}(::Type{Ptr{T}}, ::Type{T}, m::Base.RefValue{FV}) =
-         Ptr{T}(Base.unsafe_convert(Ptr{FV}, m))
+Base.cconvert(::Type{<:Ptr}, v::FieldVector) = Base.RefValue(v)
+Base.unsafe_convert(::Type{Ptr{T}}, m::Base.RefValue{FV}) where {N,T,FV<:FieldVector{N,T}} =
+    Ptr{T}(Base.unsafe_convert(Ptr{FV}, m))
diff --git a/src/SUnitRange.jl b/src/SUnitRange.jl
@@ -13,7 +13,7 @@ end
 end
 
 function check_sunitrange_params(L)
-    throw(TypeError(:SUnitRange, "type parameters must be `Int`s", Tuple{Int, Int, Int}, Tuple{typeof(a), typeof(b), typeof(c)}))
+    throw(TypeError(:SUnitRange, "type parameters must be `Int`", Tuple{Int,}, Tuple{typeof(L),}))
 end
 
 @pure SUnitRange(a::Int, b::Int) = SUnitRange{a, max(0, b - a + 1)}()
diff --git a/src/linalg.jl b/src/linalg.jl
@@ -103,7 +103,7 @@ end
 
     return quote
         @_inline_meta
-        @inbounds return similar_type(a, Size($Snew))(tuple($(exprs...)))
+        @inbounds return similar_type(a, promote_type(eltype(a), eltype(b)), Size($Snew))(tuple($(exprs...)))
     end
 end
 # TODO make these more efficient
@@ -119,7 +119,7 @@ end
 
 @generated function _hcat(::Size{Sa}, ::Size{Sb}, a::Union{StaticVector,StaticMatrix}, b::Union{StaticVector,StaticMatrix}) where {Sa, Sb}
     if Sa[1] != Sb[1]
-        error("Dimension mismatch")
+        throw(DimensionMismatch("Tried to hcat arrays of size $Sa and $Sb"))
     end
 
     exprs = vcat([:(a[$i]) for i = 1:prod(Sa)],
@@ -129,7 +129,7 @@ end
 
     return quote
         @_inline_meta
-        @inbounds return similar_type(a, Size($Snew))(tuple($(exprs...)))
+        @inbounds return similar_type(a, promote_type(eltype(a), eltype(b)), Size($Snew))(tuple($(exprs...)))
     end
 end
 # TODO make these more efficient
@@ -297,8 +297,8 @@ end
     end
 end
 
-# TODO same for `RowVector`?
 @inline Size(::Union{RowVector{T, SA}, Type{RowVector{T, SA}}}) where {T, SA <: StaticArray} = Size(1, Size(SA)[1])
+@inline Size(::Union{RowVector{T, CA}, Type{RowVector{T, CA}}} where CA <: ConjVector{<:Any, SA}) where {T, SA <: StaticArray} = Size(1, Size(SA)[1])
 @inline Size(::Union{Symmetric{T,SA}, Type{Symmetric{T,SA}}}) where {T,SA<:StaticArray} = Size(SA)
 @inline Size(::Union{Hermitian{T,SA}, Type{Hermitian{T,SA}}}) where {T,SA<:StaticArray} = Size(SA)
 
diff --git a/src/mapreduce.jl b/src/mapreduce.jl
@@ -126,14 +126,15 @@ end
 end
 
 @generated function _mapreducedim(f, op, ::Size{S}, a::StaticArray, ::Type{Val{D}}, v0) where {S, D}
-    N = ndims(a)
+    N = length(S)
     Snew = ([n==D ? 1 : S[n] for n = 1:N]...)
 
     exprs = Array{Expr}(Snew)
     itr = [1:n for n = Snew]
     for i ∈ Base.product(itr...)
         expr = :v0
         for k = 1:S[D]
+            ik = collect(i)
             ik[D] = k
             expr = :(op($expr, f(a[$(ik...)])))
         end
diff --git a/src/matrix_multiply.jl b/src/matrix_multiply.jl
@@ -36,11 +36,13 @@ promote_matprod{T1,T2}(::Type{T1}, ::Type{T2}) = typeof(zero(T1)*zero(T2) + zero
 @inline *(A::StaticMatrix, B::StaticMatrix) = _A_mul_B(Size(A), Size(B), A, B)
 @inline *(A::StaticVector, B::StaticMatrix) = *(reshape(A, Size(Size(A)[1], 1)), B)
 @inline *(A::StaticVector, B::RowVector{<:Any, <:StaticVector}) = _A_mul_B(Size(A), Size(B), A, B)
+@inline *(A::StaticVector, B::RowVector{<:Any, <:ConjVector{<:Any, <:StaticVector}}) = _A_mul_B(Size(A), Size(B), A, B)
 
 @inline A_mul_B!(dest::StaticVecOrMat, A::StaticMatrix, B::StaticVector) = _A_mul_B!(Size(dest), dest, Size(A), Size(B), A, B)
 @inline A_mul_B!(dest::StaticVecOrMat, A::StaticMatrix, B::StaticMatrix) = _A_mul_B!(Size(dest), dest, Size(A), Size(B), A, B)
 @inline A_mul_B!(dest::StaticVecOrMat, A::StaticVector, B::StaticMatrix) = A_mul_B!(dest, reshape(A, Size(Size(A)[1], 1)), B)
 @inline A_mul_B!(dest::StaticVecOrMat, A::StaticVector, B::RowVector{<:Any, <:StaticVector}) = _A_mul_B!(Size(dest), dest, Size(A), Size(B), A, B)
+@inline A_mul_B!(dest::StaticVecOrMat, A::StaticVector, B::RowVector{<:Any, <:ConjVector{<:Any, <:StaticVector}}) = _A_mul_B!(Size(dest), dest, Size(A), Size(B), A, B)
 
 #@inline *{TA<:Base.LinAlg.BlasFloat,Tb}(A::StaticMatrix{TA}, b::StaticVector{Tb})
 
@@ -93,6 +95,18 @@ end
     end
 end
 
+# complex outer product
+@generated function _A_mul_B(::Size{sa}, ::Size{sb}, a::StaticVector{<: Any, Ta}, b::RowVector{Tb, <:ConjVector{<:Any, <: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)
+        @inbounds return similar_type(b, T, Size($newsize))(tuple($(exprs...)))
+    end
+end
+
 @generated function _A_mul_B(Sa::Size{sa}, Sb::Size{sb}, a::StaticMatrix{<:Any, <:Any, Ta}, b::StaticMatrix{<:Any, <:Any, Tb}) where {sa, sb, Ta, Tb}
     # Heuristic choice for amount of codegen
     if sa[1]*sa[2]*sb[2] <= 8*8*8
diff --git a/test/MMatrix.jl b/test/MMatrix.jl
@@ -54,6 +54,8 @@
         @test ((@MMatrix eye(Float32, 2, 2))::MMatrix{2, 2, Float32}).data === (1.0f0, 0.0f0, 0.0f0, 1.0f0)
         @test isa(@MMatrix(rand(Float32, 2, 2)), MMatrix{2, 2, Float32})
         @test isa(@MMatrix(randn(Float32, 2, 2)), MMatrix{2, 2, Float32})
+
+        @test MMatrix(SMatrix{1,1,Int,1}((1,))).data == (1,)
     end
 
     @testset "Methods" begin
@@ -87,5 +89,12 @@
         m[3] = 13
         m[4] = 14
         @test m.data === (11, 12, 13, 14)
+        
+        m = @MMatrix [0 0; 0 0]
+        m[1] = Int8(11)
+        m[2] = Int8(12)
+        m[3] = Int8(13)
+        m[4] = Int8(14)
+        @test m.data === (11, 12, 13, 14)
     end
 end
diff --git a/test/MVector.jl b/test/MVector.jl
@@ -62,5 +62,12 @@
         v[2] = 12
         v[3] = 13
         @test v.data === (11, 12, 13)
+        
+        v = @MVector [1.,2.,3.]
+        v[1] = Float16(11)
+        @test v.data === (11., 2., 3.)
+
+        @test_throws BoundsError setindex!(v, 4., -1)
+        @test_throws BoundsError setindex!(v, 4., 4)
     end
 end
diff --git a/test/SUnitRange.jl b/test/SUnitRange.jl
@@ -3,4 +3,8 @@
     @test length(StaticArrays.SUnitRange(1,-10)) === 0
 
     @test StaticArrays.SUnitRange(2,4)[2] === 3
+    
+    @test_throws BoundsError StaticArrays.SUnitRange(1,3)[4]
+    @test_throws Exception StaticArrays.SUnitRange{1, -1}()
+    @test_throws TypeError StaticArrays.SUnitRange{1, 1.5}()
 end
diff --git a/test/Scalar.jl b/test/Scalar.jl
@@ -2,4 +2,7 @@
     @test Scalar(2) .* [1, 2, 3] == [2, 4, 6]
     @test Scalar([1 2; 3 4]) .+ [[1 1; 1 1], [2 2; 2 2]] == [[2 3; 4 5], [3 4; 5 6]]
     @test (Scalar(1) + Scalar(1.0))::Scalar{Float64} ≈ Scalar(2.0)
+    @test_throws ErrorException Scalar(2)[2]
+    @test Scalar(2)[] == 2
+    @test Tuple(Scalar(2)) == (2,)
 end
diff --git a/test/SizedArray.jl b/test/SizedArray.jl
@@ -1,5 +1,6 @@
 @testset "SizedArray" begin
     @testset "Inner Constructors" begin
+        @test SizedArray{Tuple{2}, Int, 1}((3, 4)).data == [3, 4]
         @test SizedArray{Tuple{2}, Int, 1}([3, 4]).data == [3, 4]
         @test SizedArray{Tuple{2, 2}, Int, 2}(collect(3:6)).data == collect(3:6)
         @test size(SizedArray{Tuple{4, 5}, Int, 2}().data) == (4, 5)
@@ -11,5 +12,19 @@
         # Bad parameters
         @test_throws Exception SizedArray{Tuple{1},Int,2}()
         @test_throws Exception SArray{Tuple{3, 4},Int,1}()
+        
+        # Parameter/input size mismatch
+        @test_throws Exception SizedArray{Tuple{1},Int,2}([2; 3])
+        @test_throws Exception SizedArray{Tuple{1},Int,2}((2, 3))
+    end
+
+    # setindex
+    sa = SizedArray{Tuple{2}, Int, 1}([3, 4])
+    sa[1] = 2
+    @test sa.data == [2, 4]
+
+    @testset "back to Array" begin
+        @test Array(SizedArray{Tuple{2}, Int, 1}([3, 4])) == [3, 4]
+        @test Vector(SizedArray{Tuple{4}, Int, 1}(collect(3:6))) == collect(3:6)
     end
 end
diff --git a/test/abstractarray.jl b/test/abstractarray.jl
@@ -55,7 +55,9 @@
         @test isa(@inferred(similar(sv, Float64, Size(3,3))), MMatrix{3,3,Float64,9})
 
         @test isa(@inferred(similar(sa)), MArray{Tuple{1,1,1},Int,3,1})
+        @test isa(@inferred(similar(SArray{Tuple{1,1,1},Int,3,1})), MArray{Tuple{1,1,1},Int,3,1})
         @test isa(@inferred(similar(sa, Float64)), MArray{Tuple{1,1,1},Float64,3,1})
+        @test isa(@inferred(similar(SArray{Tuple{1,1,1},Int,3,1}, Float64)), MArray{Tuple{1,1,1},Float64,3,1})
         @test isa(@inferred(similar(sv, Size(3,3,3))), MArray{Tuple{3,3,3},Int,3,27})
         @test isa(@inferred(similar(sv, Float64, Size(3,3,3))), MArray{Tuple{3,3,3},Float64,3,27})
     end
@@ -65,5 +67,20 @@
         @test @inferred(reshape([1,2,3,4], Size(2,2)))::SizedArray{Tuple{2,2},Int,2,1} == [1 3; 2 4]
 
         @test @inferred(vec(SMatrix{2, 2}([1 2; 3 4])))::SVector{4,Int} == [1, 3, 2, 4]
+        @test @inferred(copy(SMatrix{2, 2}([1 2; 3 4]))) == [1 2; 3 4]
+    end
+    
+    @testset "full" begin
+        m_a = randn(2,2)
+        m_a = m_a*m_a.'
+        m = SMatrix{2,2}(m_a)
+        @test @inferred(full(Symmetric(m))) == m_a
+        @test @inferred(full(Symmetric(m, :L))) == m_a
+        
+        m_a = randn(2,2) + im*randn(2,2)
+        m_a = m_a*m_a'
+        m = SMatrix{2,2}(m_a)
+        @test @inferred(full(Hermitian(m))) == m_a
+        @test @inferred(full(Hermitian(m, :L))) == m_a
     end
 end
diff --git a/test/chol.jl b/test/chol.jl
@@ -0,0 +1,25 @@
+@testset "Cholesky decomposition" begin
+    @testset "1×1" begin
+        m = @SMatrix [4.0]
+        (c,) = chol(m)
+        @test c === 2.0
+    end
+
+    @testset "2×2" for i = 1:100
+        m_a = randn(2,2)
+        #non hermitian
+        @test_throws ArgumentError chol(SMatrix{2,2}(m_a))
+        m_a = m_a*m_a'
+        m = SMatrix{2,2}(m_a)
+        @test chol(m) ≈ chol(m_a)
+    end
+
+    @testset "3×3" for i = 1:100
+        m_a = randn(3,3)
+        #non hermitian
+        @test_throws ArgumentError chol(SMatrix{3,3}(m_a))
+        m_a = m_a*m_a'
+        m = SMatrix{3,3}(m_a)
+        @test chol(m) ≈ chol(m_a)
+    end
+end
diff --git a/test/det.jl b/test/det.jl
@@ -4,8 +4,12 @@
     @test det(@SMatrix [0 1 0; 1 0 0; 0 0 1]) === -1
     m = randn(Float64, 4,4)
     @test det(SMatrix{4,4}(m)) ≈ det(m)
+    
+    @test det(@SMatrix [1 2; 0 3]) === 3
 
     # Unsigned specializations
     @test det(@SMatrix [0x00 0x01; 0x01 0x00])::Int8 == -1
     @test det(@SMatrix [0x00 0x01 0x00; 0x01 0x00 0x00; 0x00 0x00 0x01])::Int8 == -1
+
+    @test_throws DimensionMismatch det(@SMatrix [0; 1])
 end
diff --git a/test/eigen.jl b/test/eigen.jl
@@ -13,6 +13,7 @@
 
     @testset "2×2" for i = 1:100
         m_a = randn(2,2)
+        @test_throws ErrorException eigvals(SMatrix{2,2}(m_a))
         m_a = m_a*m_a'
         m = SMatrix{2,2}(m_a)
 
@@ -26,6 +27,20 @@
         @test vals::SVector ≈ vals_a
         @test eigvals(m) ≈ vals
         @test (vecs*diagm(vals)*vecs')::SMatrix ≈ m
+        
+        (vals, vecs) = eig(Hermitian(m))
+        @test vals::SVector ≈ vals_a
+        @test eigvals(Hermitian(m)) ≈ vals
+        @test eigvals(Hermitian(m, :L)) ≈ vals
+        @test (vecs*diagm(vals)*vecs')::SMatrix ≈ m
+        
+        m_d = randn(SVector{2}); m = diagm(m_d)
+        (vals, vecs) = eig(Hermitian(m))
+        @test vals::SVector ≈ sort(m_d)
+        (vals, vecs) = eig(Hermitian(m, :L))
+        @test vals::SVector ≈ sort(m_d)
+        @test eigvals(m) ≈ sort(m_d)
+        @test eigvals(Hermitian(m)) ≈ sort(m_d)
     end
 
     @testset "3×3" for i = 1:100
@@ -42,7 +57,20 @@
         (vals, vecs) = eig(Symmetric(m))
         @test vals::SVector ≈ vals_a
         @test eigvals(m) ≈ vals
+        @test eigvals(Hermitian(m)) ≈ vals
+        @test eigvals(Hermitian(m, :L)) ≈ vals
         @test (vecs*diagm(vals)*vecs')::SMatrix ≈ m
+       
+        (vals, vecs) = eig(Symmetric(m, :L))
+        @test vals::SVector ≈ vals_a
+        
+        m_d = randn(SVector{3}); m = diagm(m_d)
+        (vals, vecs) = eig(Hermitian(m))
+        @test vals::SVector ≈ sort(m_d)
+        (vals, vecs) = eig(Hermitian(m, :L))
+        @test vals::SVector ≈ sort(m_d)
+        @test eigvals(m) ≈ sort(m_d)
+        @test eigvals(Hermitian(m)) ≈ sort(m_d)
     end
 
     @testset "3x3 degenerate cases" begin
diff --git a/test/fixed_size_arrays.jl b/test/fixed_size_arrays.jl
@@ -36,6 +36,7 @@ StaticArrays.similar_type{SV <: RGB, T}(::Type{SV}, ::Type{T}, ::Size{(3,)}) = R
             (RGB{Float32}(NaN, NaN, NaN), true),
         )
         @test any(isnan, p) == r
+        @test isnan(p) == r
     end
 end
 
diff --git a/test/inv.jl b/test/inv.jl
@@ -7,4 +7,9 @@
 
     @test inv(@SMatrix [0x01 0x02; 0x02 0x01])::SMatrix ≈ [-1/3 2/3; 2/3 -1/3]
     @test inv(@SMatrix [0x01 0x02 0x00; 0x02 0x01 0x00; 0x00 0x00 0x01])::SMatrix ≈ [-1/3 2/3 0; 2/3 -1/3 0; 0 0 1]
+    
+    m = triu(randn(Float64, 4,4) + eye(4))
+    @test inv(SMatrix{4,4}(m))::StaticMatrix ≈ inv(m)
+    m = tril(randn(Float64, 4,4) + eye(4))
+    @test inv(SMatrix{4,4}(m))::StaticMatrix ≈ inv(m)
 end
diff --git a/test/linalg.jl b/test/linalg.jl
@@ -28,6 +28,7 @@
         @test @inferred(I + @SMatrix([0 1; 2 3])) === @SMatrix [1 1; 2 4]
         @test @inferred(@SMatrix([0 1; 2 3]) - I) === @SMatrix [-1 1; 2 2]
         @test @inferred(I - @SMatrix([0 1; 2 3])) === @SMatrix [1 -1; -2 -2]
+        @test_throws DimensionMismatch I + @SMatrix [0 1 4; 2 3 5]
 
         @test @inferred(@SMatrix([0 1; 2 3]) * I) === @SMatrix [0 1; 2 3]
         @test @inferred(I * @SMatrix([0 1; 2 3])) === @SMatrix [0 1; 2 3]
@@ -39,15 +40,19 @@
         @test @inferred(diagm(SVector(1,2))) === @SMatrix [1 0; 0 2]
     end
 
-    @testset "one()" begin
+    @testset "one() and zero()" begin
         @test @inferred(one(SMatrix{2,2,Int})) === @SMatrix [1 0; 0 1]
         @test @inferred(one(SMatrix{2,2})) === @SMatrix [1.0 0.0; 0.0 1.0]
         @test @inferred(one(SMatrix{2})) === @SMatrix [1.0 0.0; 0.0 1.0]
         @test @inferred(one(one(SMatrix{2,2,Int}))) === @SMatrix [1 0; 0 1]
+        @test @inferred(zero(SMatrix{2,2,Int})) === @SMatrix [0 0; 0 0]
+        @test @inferred(zero(one(SMatrix{2,2,Int}))) === @SMatrix [0 0; 0 0]
 
         @test @inferred(one(MMatrix{2,2,Int}))::MMatrix == @MMatrix [1 0; 0 1]
         @test @inferred(one(MMatrix{2,2}))::MMatrix == @MMatrix [1.0 0.0; 0.0 1.0]
         @test @inferred(one(MMatrix{2}))::MMatrix == @MMatrix [1.0 0.0; 0.0 1.0]
+
+        @test_throws ErrorException one(MMatrix{2,4})
     end
 
     @testset "eye()" begin
@@ -89,6 +94,8 @@
 
         @test @inferred(vcat(SVector(1,2,3), SVector(4,5,6))) === SVector(1,2,3,4,5,6)
         @test @inferred(hcat(SVector(1,2,3), SVector(4,5,6))) === @SMatrix [1 4; 2 5; 3 6]
+        @test_throws DimensionMismatch vcat(SVector(1,2,3), @SMatrix [1 4; 2 5])
+        @test_throws DimensionMismatch hcat(SVector(1,2,3), SVector(4,5))
 
         @test @inferred(vcat(@SMatrix([1;2;3]), SVector(4,5,6))) === @SMatrix([1;2;3;4;5;6])
         @test @inferred(vcat(SVector(1,2,3), @SMatrix([4;5;6]))) === @SMatrix([1;2;3;4;5;6])
@@ -100,6 +107,9 @@
 
         @test @inferred(vcat(SVector(1),SVector(2),SVector(3),SVector(4))) === SVector(1,2,3,4)
         @test @inferred(hcat(SVector(1),SVector(2),SVector(3),SVector(4))) === SMatrix{1,4}(1,2,3,4)
+
+        vcat(SVector(1.0f0), SVector(1.0)) === SVector(1.0, 1.0)
+        hcat(SVector(1.0f0), SVector(1.0)) === SMatrix{1,2}(1.0, 1.0)
     end
 
     @testset "normalization" begin
@@ -113,9 +123,11 @@
         @test vecnorm(@SMatrix [1 2; 3 4.0+im]) ≈ vecnorm([1 2; 3 4.0+im])
 
         @test normalize(SVector(1,2,3)) ≈ normalize([1,2,3])
+        @test normalize(SVector(1,2,3), 1) ≈ normalize([1,2,3], 1)
     end
 
     @testset "trace" begin
         @test trace(@SMatrix [1.0 2.0; 3.0 4.0]) === 5.0
+        @test_throws DimensionMismatch trace(@SMatrix rand(5,4))
     end
 end
diff --git a/test/mapreduce.jl b/test/mapreduce.jl
diff --git a/test/matrix_multiply.jl b/test/matrix_multiply.jl
diff --git a/test/runtests.jl b/test/runtests.jl
