From 69da0632e77e764b1660b4b920f8d3ef04823ec5 Mon Sep 17 00:00:00 2001 From: Tianyi Pu <912396513@qq.com> Date: Sat, 25 Mar 2023 23:03:43 +0000 Subject: [PATCH 01/11] add unary and number broadcast; add `_all` --- src/ToeplitzMatrices.jl | 7 ++++--- src/special.jl | 10 ++++------ src/toeplitz.jl | 7 ++++--- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/ToeplitzMatrices.jl b/src/ToeplitzMatrices.jl index 687cd4f..b4abf44 100644 --- a/src/ToeplitzMatrices.jl +++ b/src/ToeplitzMatrices.jl @@ -2,8 +2,9 @@ module ToeplitzMatrices # import StatsBase: levinson!, levinson 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: adjoint, convert, transpose, size, getindex, similar, copy, getproperty, inv, sqrt, copyto!, reverse, conj, zero, fill!, checkbounds, real, isfinite, DimsInteger, _all import Base: ==, +, -, *, \ +import Base.Broadcast: broadcasted, DefaultMatrixStyle import LinearAlgebra: Cholesky, Factorization import LinearAlgebra: ldiv!, factorize, lmul!, pinv, eigvals, cholesky!, cholesky, tril!, triu!, checksquare, rmul!, dot, mul!, tril, triu import LinearAlgebra: UpperTriangular, LowerTriangular, Symmetric, Adjoint @@ -16,7 +17,7 @@ export durbin, trench, levinson 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 @@ -28,7 +29,7 @@ AbstractArray{T}(A::AbstractToeplitz) where T = AbstractToeplitz{T}(A) convert(::Type{AbstractToeplitz{T}}, A::AbstractToeplitz) where T = AbstractToeplitz{T}(A) isconcrete(A::AbstractToeplitz) = isconcretetype(typeof(A.vc)) && isconcretetype(typeof(A.vr)) -iszero(A::AbstractToeplitz) = iszero(A.vc) && iszero(A.vr) +_all(f, A::AbstractToeplitz, ::Colon) = _all(f, A.vc, :) && _all(f, A.vr, :) """ ToeplitzFactorization diff --git a/src/special.jl b/src/special.jl index 6e271c8..e78c2e7 100644 --- a/src/special.jl +++ b/src/special.jl @@ -14,9 +14,10 @@ for TYPE in (:SymmetricToeplitz, :Circulant, :LowerTriangularToeplitz, :UpperTri size(A::$TYPE) = (length(A.v),length(A.v)) + broadcasted(::DefaultMatrixStyle, f, A::$TYPE) = $TYPE(f.(A.v)) + broadcasted(::DefaultMatrixStyle, f, x::Number, A::$TYPE) = $TYPE(f.(x, A.v)) + broadcasted(::DefaultMatrixStyle, f, A::$TYPE, x::Number) = $TYPE(f.(A.v, x)) adjoint(A::$TYPE) = transpose(conj(A)) - (*)(scalar::Number, C::$TYPE) = $TYPE(scalar * C.v) - (*)(C::$TYPE, scalar::Number) = $TYPE(C.v * scalar) (==)(A::$TYPE,B::$TYPE) = A.v==B.v function zero!(A::$TYPE) if isconcrete(A) @@ -40,12 +41,9 @@ for TYPE in (:SymmetricToeplitz, :Circulant, :LowerTriangularToeplitz, :UpperTri A end end - for fun in (:zero, :conj, :copy, :-, :real, :imag) + for fun in (:zero, :copy) @eval $fun(A::$TYPE) = $TYPE($fun(A.v)) end - for fun in (:iszero,) - @eval $fun(A::$TYPE) = $fun(A.v) - end for op in (:+, :-) @eval $op(A::$TYPE,B::$TYPE) = $TYPE($op(A.v,B.v)) end diff --git a/src/toeplitz.jl b/src/toeplitz.jl index 5bb97fd..5eb2cd8 100644 --- a/src/toeplitz.jl +++ b/src/toeplitz.jl @@ -53,6 +53,9 @@ Base.@propagate_inbounds function getindex(A::AbstractToeplitz, i::Integer, j::I end end +broadcasted(::DefaultMatrixStyle, f, A::AbstractToeplitz) = Toeplitz(f.(A.vc), f.(A.vr)) +broadcasted(::DefaultMatrixStyle, f, x::Number, A::AbstractToeplitz) = Toeplitz(f.(x, A.vc), f.(x, A.vr)) + 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 +113,7 @@ function similar(A::AbstractToeplitz, T::Type, dims::Dims{2}) vr[1]=vc[1] 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 +132,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) From 5a619d5f86964b95246ff6f9fe7ac4e20e27a8f0 Mon Sep 17 00:00:00 2001 From: Tianyi Pu <912396513@qq.com> Date: Sun, 26 Mar 2023 11:22:19 +0100 Subject: [PATCH 02/11] Update toeplitz.jl --- src/toeplitz.jl | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/src/toeplitz.jl b/src/toeplitz.jl index 5eb2cd8..5536e2b 100644 --- a/src/toeplitz.jl +++ b/src/toeplitz.jl @@ -9,7 +9,7 @@ struct Toeplitz{T, VC<:AbstractVector{T}, VR<:AbstractVector{T}} <: AbstractToep vr::VR function Toeplitz{T, VC, VR}(vc::VC, vr::VR) where {T, VC<:AbstractVector{T}, VR<:AbstractVector{T}} - 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) @@ -42,6 +42,22 @@ AbstractToeplitz{T}(A::Toeplitz) where T = Toeplitz{T}(A) convert(::Type{Toeplitz{T}}, A::AbstractToeplitz) where {T} = Toeplitz{T}(A) convert(::Type{Toeplitz}, A::AbstractToeplitz) = Toeplitz(A) +# Dims for disambiguity +for DIM in (:DimsInteger, :Dims) + @eval function similar(::Type{Toeplitz{T, VC, VR}}, dims::$DIM{2}) where {T, VC, VR} + vc = similar(VC, dims[1]) + vr = similar(VR, dims[2]) + vr[1] = vc[1] + Toeplitz(vc, vr) + end + @eval function similar(A::AbstractToeplitz, T::Type = eltype(A), dims::$DIM{2} = size(A)) + vc = similar(A.vc, T, dims[1]) + vr = similar(A.vr, T, dims[2]) + vr[1] = vc[1] + Toeplitz{T}(vc, vr) + end +end + # Retrieve an entry Base.@propagate_inbounds function getindex(A::AbstractToeplitz, i::Integer, j::Integer) @boundscheck checkbounds(A,i,j) @@ -107,12 +123,6 @@ end adjoint(A::AbstractToeplitz) = transpose(conj(A)) transpose(A::AbstractToeplitz) = Toeplitz(A.vr, A.vc) -function similar(A::AbstractToeplitz, T::Type, dims::Dims{2}) - vc=similar(A.vc, T, dims[1]) - vr=similar(A.vr, T, dims[2]) - vr[1]=vc[1] - Toeplitz{T}(vc,vr) -end for fun in (:zero, :copy) @eval $fun(A::AbstractToeplitz)=Toeplitz($fun(A.vc),$fun(A.vr)) end From 34e56a56c40a9dc09fb8457516805bd7747ee21b Mon Sep 17 00:00:00 2001 From: Tianyi Pu <912396513@qq.com> Date: Sun, 26 Mar 2023 11:30:49 +0100 Subject: [PATCH 03/11] Update toeplitz.jl --- src/toeplitz.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/toeplitz.jl b/src/toeplitz.jl index 5536e2b..ffc7e46 100644 --- a/src/toeplitz.jl +++ b/src/toeplitz.jl @@ -71,6 +71,7 @@ end broadcasted(::DefaultMatrixStyle, f, A::AbstractToeplitz) = Toeplitz(f.(A.vc), f.(A.vr)) broadcasted(::DefaultMatrixStyle, f, x::Number, A::AbstractToeplitz) = Toeplitz(f.(x, A.vc), f.(x, A.vr)) +broadcasted(::DefaultMatrixStyle, 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")) From 797625c12b9ae1d4cf8e014059d3289230bdfe67 Mon Sep 17 00:00:00 2001 From: Tianyi Pu <912396513@qq.com> Date: Sun, 26 Mar 2023 11:42:50 +0100 Subject: [PATCH 04/11] Update runtests.jl --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index df1b49f..bfebdef 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -362,6 +362,7 @@ end T=copy(TA) end @test fill!(Toeplitz(zeros(2,2)),1) == ones(2,2) + @test isa(similar(Toeplitz{Int, Vector{Int}, Vector{Int}}, 5, 3), Toeplitz) @testset "aliasing" begin v = [1,2,3] From e9e57f38e3ea2decf3c2e5ba398990387779dd69 Mon Sep 17 00:00:00 2001 From: Tianyi Pu <912396513@qq.com> Date: Sun, 26 Mar 2023 12:12:09 +0100 Subject: [PATCH 05/11] add a case where broadcast breaks --- test/runtests.jl | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index bfebdef..e68eb38 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -378,6 +378,17 @@ 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) + + @test U .+ 1 == Matrix(U) .+ 1 +end + @testset "Circulant mathematics" begin C1 = Circulant(rand(5)) C2 = Circulant(rand(5)) From 51900b6f9a92980f46e2bde2dc68f9371e9cd1a8 Mon Sep 17 00:00:00 2001 From: Tianyi Pu <912396513@qq.com> Date: Sun, 26 Mar 2023 13:44:12 +0100 Subject: [PATCH 06/11] fixes and more tests --- README.md | 11 ++++------- src/special.jl | 19 ++++++++++++++++++- src/toeplitz.jl | 9 ++++++--- test/runtests.jl | 8 +++++++- 4 files changed, 35 insertions(+), 12 deletions(-) diff --git a/README.md b/README.md index 8915c56..b2b9ce3 100644 --- a/README.md +++ b/README.md @@ -78,10 +78,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!|✓|✗|✗|✓|✓|✗| @@ -90,19 +87,19 @@ The `reverse` can transform between Hankel and Toeplitz. It is used to achieve f |triu|✓|✓|✓|✓|✓|✗| |+|✓|✓|✓|✓|✓|✓| |-|✓|✓|✓|✓|✓|✓| -|scalar
mult|✓|✓|✓|✓|✓|✓| |==|✓|✓|✓|✓|✓|✓| |issymmetric||||||| |istriu||||||| |istril||||||| -|iszero|✓|✓|✓|✓|✓|| |isone||||||| |copyto!|✓|✓|✓|✓|✓|✓| |reverse|✓|✓|✓|✓|✓|✓| -|broadcast||||||| +|\_all|✓|✓|✓|✓|✓|| +|unary broadcast|✓|✓|✓|✓|✓|| +|number broadcast|✓|✓|✓|✓|✓|| |broadcast!||||||| -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`. diff --git a/src/special.jl b/src/special.jl index e78c2e7..a8a31c2 100644 --- a/src/special.jl +++ b/src/special.jl @@ -51,7 +51,8 @@ for TYPE in (:SymmetricToeplitz, :Circulant, :LowerTriangularToeplitz, :UpperTri @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) @@ -102,6 +103,22 @@ transpose(A::Circulant) = Circulant(A.vr) 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) diff --git a/src/toeplitz.jl b/src/toeplitz.jl index ffc7e46..d86938a 100644 --- a/src/toeplitz.jl +++ b/src/toeplitz.jl @@ -69,9 +69,12 @@ Base.@propagate_inbounds function getindex(A::AbstractToeplitz, i::Integer, j::I end end -broadcasted(::DefaultMatrixStyle, f, A::AbstractToeplitz) = Toeplitz(f.(A.vc), f.(A.vr)) -broadcasted(::DefaultMatrixStyle, f, x::Number, A::AbstractToeplitz) = Toeplitz(f.(x, A.vc), f.(x, A.vr)) -broadcasted(::DefaultMatrixStyle, f, A::AbstractToeplitz, x::Number) = Toeplitz(f.(A.vc, x), f.(A.vr, x)) +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")) diff --git a/test/runtests.jl b/test/runtests.jl index e68eb38..ad1c99f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -386,7 +386,13 @@ end U = UpperTriangular(T) L = LowerTriangular(T) - @test U .+ 1 == Matrix(U) .+ 1 + 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 From 652caaef8deb9b5443277f7e6121e3ca9fbf282c Mon Sep 17 00:00:00 2001 From: Tianyi Pu <912396513@qq.com> Date: Sun, 26 Mar 2023 13:58:29 +0100 Subject: [PATCH 07/11] fix --- src/special.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/special.jl b/src/special.jl index a8a31c2..15f323e 100644 --- a/src/special.jl +++ b/src/special.jl @@ -14,9 +14,6 @@ for TYPE in (:SymmetricToeplitz, :Circulant, :LowerTriangularToeplitz, :UpperTri size(A::$TYPE) = (length(A.v),length(A.v)) - broadcasted(::DefaultMatrixStyle, f, A::$TYPE) = $TYPE(f.(A.v)) - broadcasted(::DefaultMatrixStyle, f, x::Number, A::$TYPE) = $TYPE(f.(x, A.v)) - broadcasted(::DefaultMatrixStyle, f, A::$TYPE, x::Number) = $TYPE(f.(A.v, x)) adjoint(A::$TYPE) = transpose(conj(A)) (==)(A::$TYPE,B::$TYPE) = A.v==B.v function zero!(A::$TYPE) From 3f3fb59e0237992a20cb050ddf231c82fe7b79a2 Mon Sep 17 00:00:00 2001 From: Tianyi Pu <912396513@qq.com> Date: Mon, 27 Mar 2023 17:54:54 +0100 Subject: [PATCH 08/11] hankel --- README.md | 6 +++--- src/hankel.jl | 8 +++++--- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index b2b9ce3..f0212cb 100644 --- a/README.md +++ b/README.md @@ -94,9 +94,9 @@ The `reverse` can transform between Hankel and Toeplitz. It is used to achieve f |isone||||||| |copyto!|✓|✓|✓|✓|✓|✓| |reverse|✓|✓|✓|✓|✓|✓| -|\_all|✓|✓|✓|✓|✓|| -|unary broadcast|✓|✓|✓|✓|✓|| -|number broadcast|✓|✓|✓|✓|✓|| +|\_all|✓|✓|✓|✓|✓|✓| +|unary broadcast|✓|✓|✓|✓|✓|✓| +|number broadcast|✓|✓|✓|✓|✓|✓| |broadcast!||||||| Note that `+` and `-` could be removed once binary `broadcast` is implemented. diff --git a/src/hankel.jl b/src/hankel.jl index 6d12ee4..dd8685c 100644 --- a/src/hankel.jl +++ b/src/hankel.jl @@ -59,6 +59,10 @@ end 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 @@ -70,7 +74,7 @@ Base.@propagate_inbounds function getindex(A::Hankel, i::Integer, j::Integer) end similar(A::Hankel, T::Type, dims::DimsInteger{2}) = Hankel{T}(similar(A.v, T, dims[1]+dims[2]-true), dims) similar(A::Hankel, T::Type, dims::Tuple{Int64,Int64}) = Hankel{T}(similar(A.v, T, dims[1]+dims[2]-true), dims) # for ambiguity with `similar(a::AbstractArray, ::Type{T}, dims::Tuple{Vararg{Int64, N}}) where {T, N}` in Base -for fun in (:zero, :conj, :copy, :-, :similar, :real, :imag) +for fun in (:zero, :copy, :similar) @eval $fun(A::Hankel) = Hankel($fun(A.v), size(A)) end for op in (:+, :-) @@ -100,8 +104,6 @@ end 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) From 5da10122f1186463e376200acbe5eb7231d2271d Mon Sep 17 00:00:00 2001 From: Tianyi Pu <912396513@qq.com> Date: Mon, 27 Mar 2023 18:19:16 +0100 Subject: [PATCH 09/11] Update hankel.jl --- src/hankel.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/hankel.jl b/src/hankel.jl index dd8685c..42db9ed 100644 --- a/src/hankel.jl +++ b/src/hankel.jl @@ -72,9 +72,9 @@ Base.@propagate_inbounds function getindex(A::Hankel, i::Integer, j::Integer) @boundscheck checkbounds(A, i, j) return A.v[i+j-1] end -similar(A::Hankel, T::Type, dims::DimsInteger{2}) = Hankel{T}(similar(A.v, T, dims[1]+dims[2]-true), dims) -similar(A::Hankel, T::Type, dims::Tuple{Int64,Int64}) = Hankel{T}(similar(A.v, T, dims[1]+dims[2]-true), dims) # for ambiguity with `similar(a::AbstractArray, ::Type{T}, dims::Tuple{Vararg{Int64, N}}) where {T, N}` in Base -for fun in (:zero, :copy, :similar) +similar(A::Hankel, T::Type = eltype(A), dims::DimsInteger{2} = size(A)) = Hankel{T}(similar(A.v, T, dims[1]+dims[2]-true), dims) +similar(A::Hankel, T::Type = eltype(A), dims::Tuple{Int64,Int64} = size(A)) = Hankel{T}(similar(A.v, T, dims[1]+dims[2]-true), dims) # for ambiguity with `similar(a::AbstractArray, ::Type{T}, dims::Tuple{Vararg{Int64, N}}) where {T, N}` in Base +for fun in (:zero, :copy) @eval $fun(A::Hankel) = Hankel($fun(A.v), size(A)) end for op in (:+, :-) From 343ace00625be75b0b0df2b813cbdd628c0d875d Mon Sep 17 00:00:00 2001 From: Tianyi Pu <912396513@qq.com> Date: Mon, 27 Mar 2023 22:11:33 +0100 Subject: [PATCH 10/11] `_all` to `all` --- src/ToeplitzMatrices.jl | 4 ++-- src/hankel.jl | 2 +- src/special.jl | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/ToeplitzMatrices.jl b/src/ToeplitzMatrices.jl index b4abf44..d939dd8 100644 --- a/src/ToeplitzMatrices.jl +++ b/src/ToeplitzMatrices.jl @@ -2,7 +2,7 @@ module ToeplitzMatrices # import StatsBase: levinson!, levinson import DSP: conv -import Base: adjoint, convert, transpose, size, getindex, similar, copy, getproperty, inv, sqrt, copyto!, reverse, conj, zero, fill!, checkbounds, real, isfinite, DimsInteger, _all +import Base: adjoint, convert, transpose, size, getindex, similar, copy, getproperty, inv, sqrt, copyto!, reverse, conj, zero, fill!, checkbounds, real, isfinite, DimsInteger, all import Base: ==, +, -, *, \ import Base.Broadcast: broadcasted, DefaultMatrixStyle import LinearAlgebra: Cholesky, Factorization @@ -29,7 +29,7 @@ AbstractArray{T}(A::AbstractToeplitz) where T = AbstractToeplitz{T}(A) convert(::Type{AbstractToeplitz{T}}, A::AbstractToeplitz) where T = AbstractToeplitz{T}(A) isconcrete(A::AbstractToeplitz) = isconcretetype(typeof(A.vc)) && isconcretetype(typeof(A.vr)) -_all(f, A::AbstractToeplitz, ::Colon) = _all(f, A.vc, :) && _all(f, A.vr, :) +all(f, A::AbstractToeplitz, ::Colon) = all(f, A.vc, :) && all(f, A.vr, :) """ ToeplitzFactorization diff --git a/src/hankel.jl b/src/hankel.jl index 42db9ed..5a2b4fb 100644 --- a/src/hankel.jl +++ b/src/hankel.jl @@ -62,7 +62,7 @@ convert(::Type{Hankel{T}}, A::Hankel) where {T<:Number} = Hankel{T}(convert(Abst 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, :) +all(f, A::Hankel, ::Colon) = all(f, A.v, :) # Size size(H::Hankel) = H.size diff --git a/src/special.jl b/src/special.jl index 15f323e..2948b5e 100644 --- a/src/special.jl +++ b/src/special.jl @@ -101,8 +101,8 @@ 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, :) +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) From 4ed4ab0fd3b41dc24076d00623af0dab87568d25 Mon Sep 17 00:00:00 2001 From: Tianyi Pu <912396513@qq.com> Date: Tue, 14 Nov 2023 10:44:35 +0000 Subject: [PATCH 11/11] resolve conflict and add mapreduce for + and max --- src/ToeplitzMatrices.jl | 7 ++++++- src/toeplitz.jl | 16 ---------------- test/runtests.jl | 1 - 3 files changed, 6 insertions(+), 18 deletions(-) diff --git a/src/ToeplitzMatrices.jl b/src/ToeplitzMatrices.jl index a30d440..fd5730b 100644 --- a/src/ToeplitzMatrices.jl +++ b/src/ToeplitzMatrices.jl @@ -1,7 +1,8 @@ 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 @@ -47,6 +48,10 @@ AbstractArray{T}(A::AbstractToeplitz) where T = AbstractToeplitz{T}(A) 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) diff --git a/src/toeplitz.jl b/src/toeplitz.jl index b89f2ea..ca3216f 100644 --- a/src/toeplitz.jl +++ b/src/toeplitz.jl @@ -43,22 +43,6 @@ AbstractToeplitz{T}(A::Toeplitz) where T = Toeplitz{T}(A) convert(::Type{Toeplitz{T}}, A::AbstractToeplitz) where {T} = Toeplitz{T}(A) convert(::Type{Toeplitz}, A::AbstractToeplitz) = Toeplitz(A) -# Dims for disambiguity -for DIM in (:DimsInteger, :Dims) - @eval function similar(::Type{Toeplitz{T, VC, VR}}, dims::$DIM{2}) where {T, VC, VR} - vc = similar(VC, dims[1]) - vr = similar(VR, dims[2]) - vr[1] = vc[1] - Toeplitz(vc, vr) - end - @eval function similar(A::AbstractToeplitz, T::Type = eltype(A), dims::$DIM{2} = size(A)) - vc = similar(A.vc, T, dims[1]) - vr = similar(A.vr, T, dims[2]) - vr[1] = vc[1] - Toeplitz{T}(vc, vr) - end -end - # Retrieve an entry Base.@propagate_inbounds function getindex(A::AbstractToeplitz, i::Integer, j::Integer) @boundscheck checkbounds(A,i,j) diff --git a/test/runtests.jl b/test/runtests.jl index 3a907dc..c1f0c55 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -377,7 +377,6 @@ end T=copy(TA) end @test fill!(Toeplitz(zeros(2,2)),1) == ones(2,2) - @test isa(similar(Toeplitz{Int, Vector{Int}, Vector{Int}}, 5, 3), Toeplitz) @testset "diag" begin H = Hankel(1:11, 4, 8)