add unary and number broadcast; add all

README.md

@@ -82,10 +82,7 @@ The `reverse` can transform between Hankel and Toeplitz. It is used to achieve f
 |copy|✓|✓|✓|✓|✓|✓|
 |similar|✓|✓|✓|✓|✓|✓|
 |zero|✓|✓|✓|✓|✓|✓|
-|real|✓|✓|✓|✓|✓|✓|
-|imag|✓|✓|✓|✓|✓|✓|
 |fill!|✓|✗|✗|✗|✗|✓|
-|conj|✓|✓|✓|✓|✓|✓|
 |transpose|✓|✓|✓|✓|✓|✓|
 |adjoint|✓|✓|✓|✓|✓|✓|
 |tril!|✓|✗|✗|✓|✓|✗|
@@ -94,22 +91,22 @@ The `reverse` can transform between Hankel and Toeplitz. It is used to achieve f
 |triu|✓|✓|✓|✓|✓|✗|
 |+|✓|✓|✓|✓|✓|✓|
 |-|✓|✓|✓|✓|✓|✓|
-|scalar<br>mult|✓|✓|✓|✓|✓|✓|
 |==|✓|✓|✓|✓|✓|✓|
 |issymmetric|||||||
 |istriu|||||||
 |istril|||||||
-|iszero|✓|✓|✓|✓|✓||
 |isone|||||||
 |diag|✓|✓|✓|✓|✓|✓|
 |copyto!|✓|✓|✓|✓|✓|✓|
 |reverse|✓|✓|✓|✓|✓|✓|
-|broadcast|||||||
+|\_all|✓|✓|✓|✓|✓|✓|
+|unary broadcast|✓|✓|✓|✓|✓|✓|
+|number broadcast|✓|✓|✓|✓|✓|✓|
 |broadcast!|||||||
 
 </details>
 
-Note that scalar multiplication, `conj`, `+` and `-` could be removed once `broadcast` is implemented.
+Note that `+` and `-` could be removed once binary `broadcast` is implemented.
 
 `reverse(Hankel)` returns a `Toeplitz`, while `reverse(AbstractToeplitz)` returns a `Hankel`.
 

src/ToeplitzMatrices.jl

@@ -1,10 +1,12 @@
 module ToeplitzMatrices
 import DSP: conv
 
-import Base: adjoint, convert, transpose, size, getindex, similar, copy, getproperty, inv, sqrt, copyto!, reverse, conj, zero, fill!, checkbounds, real, imag, isfinite, DimsInteger, iszero
+import Base: convert, size, getindex, similar, copy, getproperty, copyto!, reverse, zero, fill!, checkbounds, isfinite, DimsInteger, iszero
+import Base: all, adjoint, transpose, real, imag, inv, sqrt, conj
 import Base: parent
 import Base: ==, +, -, *, \
 import Base: AbstractMatrix
+import Base.Broadcast: broadcasted, DefaultMatrixStyle
 import LinearAlgebra: Cholesky, Factorization
 import LinearAlgebra: ldiv!, factorize, lmul!, pinv, eigvals, eigvecs, eigen, Eigen, det
 import LinearAlgebra: cholesky!, cholesky, tril!, triu!, checksquare, rmul!, dot, mul!, tril, triu, diag
@@ -34,7 +36,7 @@
 include("iterativeLinearSolvers.jl")
 
 # Abstract
-abstract type AbstractToeplitz{T<:Number} <: AbstractMatrix{T} end
+abstract type AbstractToeplitz{T} <: AbstractMatrix{T} end
 
 size(A::AbstractToeplitz) = (length(A.vc),length(A.vr))
 @inline _vr(A::AbstractToeplitz) = A.vr
@@ -46,6 +48,11 @@
 convert(::Type{AbstractToeplitz{T}}, A::AbstractToeplitz) where T = AbstractToeplitz{T}(A)
 
 isconcrete(A::AbstractToeplitz) = isconcretetype(typeof(A.vc)) && isconcretetype(typeof(A.vr))
+
+for op in (:+, :max)
+    @eval mapreduce(f, ::typeof($op), A::AbstractToeplitz) = mapreduce(identity, $op, [mapreduce(f, $op, diag(A, k)) for k in -size(A,1)+1:size(A,2)-1])
+end
+all(f, A::AbstractToeplitz, ::Colon) = all(f, A.vc, :) && all(f, A.vr, :)
 iszero(A::AbstractToeplitz) = iszero(A.vc) && iszero(A.vr)
 function diag(A::AbstractToeplitz, k::Integer=0)
     if k >= size(A, 2) || -k >= size(A, 1)

src/hankel.jl

@@ -60,6 +60,10 @@
 convert(::Type{AbstractArray{T}}, A::Hankel) where {T<:Number} = convert(Hankel{T}, A)
 convert(::Type{AbstractMatrix{T}}, A::Hankel) where {T<:Number} = convert(Hankel{T}, A)
 convert(::Type{Hankel{T}}, A::Hankel) where {T<:Number} = Hankel{T}(convert(AbstractVector{T}, A.v), size(A))
+broadcasted(::DefaultMatrixStyle, f, A::Hankel) = Hankel(f.(A.v), A.size)
+broadcasted(::DefaultMatrixStyle, f, x::Number, A::Hankel) = Hankel(f.(x, A.v), A.size)
+broadcasted(::DefaultMatrixStyle, f, A::Hankel, x::Number) = Hankel(f.(A.v, x), A.size)
+all(f, A::Hankel, ::Colon) = all(f, A.v, :)
 
 # Size
 size(H::Hankel) = H.size
@@ -72,7 +76,7 @@
     return A.v[i+j-1]
 end
 AbstractMatrix{T}(A::Hankel) where {T} = Hankel{T}(AbstractVector{T}(A.v), A.size)
-for fun in (:zero, :conj, :copy, :-, :similar, :real, :imag)
+for fun in (:zero, :copy)
     @eval $fun(A::Hankel) = Hankel($fun(A.v), size(A))
 end
 for op in (:+, :-)
@@ -102,8 +106,6 @@
 transpose(A::Hankel) = Hankel(A.v, reverse(size(A)))
 adjoint(A::Hankel) = transpose(conj(A))
 (==)(A::Hankel, B::Hankel) = A.v == B.v && size(A) == size(B)
-(*)(scalar::Number, C::Hankel) = Hankel(scalar * C.v, size(C))
-(*)(C::Hankel,scalar::Number) = Hankel(C.v * scalar, size(C))
 
 isconcrete(A::Hankel) = isconcretetype(A.v)
 

src/special.jl

@@ -74,7 +74,8 @@
         @eval $TYPE(v::$TY) = $TYPE{eltype(v)}(v)
     end
 end
-TriangularToeplitz{T}=Union{UpperTriangularToeplitz{T},LowerTriangularToeplitz{T}}
+SymToeplitzOrCirc{T} = Union{SymmetricToeplitz{T}, Circulant{T}}
+TriangularToeplitz{T} = Union{UpperTriangularToeplitz{T}, LowerTriangularToeplitz{T}}
 
 # vc and vr
 function getproperty(A::SymmetricToeplitz, s::Symbol)
@@ -129,6 +130,22 @@
 transpose(A::LowerTriangularToeplitz) = UpperTriangularToeplitz(A.v)
 transpose(A::UpperTriangularToeplitz) = LowerTriangularToeplitz(A.v)
 
+# _all
+all(f, A::SymToeplitzOrCirc, ::Colon) = all(f, A.v, :)
+all(f, A::TriangularToeplitz, ::Colon) = f(zero(eltype(A))) && all(f, A.v, :)
+
+# broadcast
+for TYPE in (:SymmetricToeplitz, :Circulant)
+    @eval broadcasted(::DefaultMatrixStyle, f, A::$TYPE) = $TYPE(f.(A.v))
+    @eval broadcasted(::DefaultMatrixStyle, f, x::Number, A::$TYPE) = $TYPE(f.(x, A.v))
+    @eval broadcasted(::DefaultMatrixStyle, f, A::$TYPE, x::Number) = $TYPE(f.(A.v, x))
+end
+for TYPE in (:UpperTriangularToeplitz, :LowerTriangularToeplitz)
+    @eval broadcasted(::DefaultMatrixStyle, f, A::$TYPE) = iszero(f(zero(eltype(A)))) ? $TYPE(f.(A.v)) : _toep_broadcast(f, A)
+    @eval broadcasted(::DefaultMatrixStyle, f, x::Number, A::$TYPE) = iszero(f(x, zero(eltype(A)))) ? $TYPE(f.(x, A.v)) : _toep_broadcast(f, x, A)
+    @eval broadcasted(::DefaultMatrixStyle, f, A::$TYPE, x::Number) = iszero(f(zero(eltype(A)), x)) ? $TYPE(f.(A.v, x)) : _toep_broadcast(f, A, x)
+end
+
 # getindex
 Base.@propagate_inbounds function getindex(A::SymmetricToeplitz, i::Integer, j::Integer)
     @boundscheck checkbounds(A,i,j)

src/toeplitz.jl

@@ -10,7 +10,7 @@ struct Toeplitz{T, VC<:AbstractVector{T}, VR<:AbstractVector{T}} <: AbstractToep
 
     function Toeplitz{T, VC, VR}(vc::VC, vr::VR) where {T, VC<:AbstractVector{T}, VR<:AbstractVector{T}}
         require_one_based_indexing(vr, vc)
-        if first(vc) != first(vr)
+        if !isequal(first(vc), first(vr))
             error("First element of the vectors must be the same")
         end
         return new{T,VC,VR}(vc, vr)
@@ -54,6 +54,13 @@ Base.@propagate_inbounds function getindex(A::AbstractToeplitz, i::Integer, j::I
     end
 end
 
+broadcasted(::DefaultMatrixStyle, f, A::AbstractToeplitz) = _toep_broadcast(f, A)
+broadcasted(::DefaultMatrixStyle, f, x::Number, A::AbstractToeplitz) = _toep_broadcast(f, x, A)
+broadcasted(::DefaultMatrixStyle, f, A::AbstractToeplitz, x::Number) = _toep_broadcast(f, A, x)
+_toep_broadcast(f, A::AbstractToeplitz) = Toeplitz(f.(A.vc), f.(A.vr))
+_toep_broadcast(f, x::Number, A::AbstractToeplitz) = Toeplitz(f.(x, A.vc), f.(x, A.vr))
+_toep_broadcast(f, A::AbstractToeplitz, x::Number) = Toeplitz(f.(A.vc, x), f.(A.vr, x))
+
 checknonaliased(A::Toeplitz) = Base.mightalias(A.vc, A.vr) && throw(ArgumentError("Cannot modify Toeplitz matrices in place with aliased data"))
 
 function tril!(A::Toeplitz, k::Integer=0)
@@ -110,7 +117,7 @@ function AbstractMatrix{T}(A::AbstractToeplitz) where {T}
     vr = AbstractVector{T}(_vr(A))
     Toeplitz{T}(vc,vr)
 end
-for fun in (:zero, :conj, :copy, :-, :real, :imag)
+for fun in (:zero, :copy)
     @eval $fun(A::AbstractToeplitz)=Toeplitz($fun(A.vc),$fun(A.vr))
 end
 for op in (:+, :-)
@@ -129,8 +136,6 @@ function fill!(A::Toeplitz, x::Number)
     fill!(A.vr,x)
     A
 end
-(*)(scalar::Number, C::AbstractToeplitz) = Toeplitz(scalar * C.vc, scalar * C.vr)
-(*)(C::AbstractToeplitz,scalar::Number) = Toeplitz(C.vc * scalar, C.vr * scalar)
 
 function lmul!(x::Number, A::Toeplitz)
     checknonaliased(A)

test/runtests.jl

@@ -430,6 +430,23 @@ end
     end
 end
 
+@testset "Broadcast" begin
+    A = rand(ComplexF64, 3, 3)
+    T = Toeplitz(A)
+    S = Symmetric(T)
+    C = Circulant(T)
+    U = UpperTriangular(T)
+    L = LowerTriangular(T)
+
+    for M in (T, S, C, U, L)
+        @testset "$(typeof(M))" begin
+            @test M .+ 1 == Matrix(M) .+ 1
+            @test exp.(M) == exp.(Matrix(M))
+            @test isa(M .* 2, typeof(M))
+        end
+    end
+end
+
 @testset "Circulant mathematics" begin
     C1 = Circulant(rand(5))
     C2 = Circulant(rand(5))