Merge branch 'master' into julia1.6

Author: jishnub

.github/workflows/CI.yml

@@ -7,7 +7,18 @@ on:
   pull_request:
 
 jobs:
+  pre_job:
+    # continue-on-error: true # Uncomment once integration is finished
+    runs-on: ubuntu-latest
+    # Map a step output to a job output
+    outputs:
+      should_skip: ${{ steps.skip_check.outputs.should_skip }}
+    steps:
+      - id: skip_check
+        uses: fkirc/skip-duplicate-actions@v5
   test:
+    needs: pre_job
+    if: needs.pre_job.outputs.should_skip != 'true'
     name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
     runs-on: ${{ matrix.os }}
     strategy:

Project.toml

@@ -1,6 +1,6 @@
 name = "ToeplitzMatrices"
 uuid = "c751599d-da0a-543b-9d20-d0a503d91d24"
-version = "0.9.0"
+version = "0.8.2"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
@@ -10,8 +10,11 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
 AbstractFFTs = "0.4, 0.5, 1"
+Aqua = "0.6"
 DSP = "0.7.7"
 StatsBase = "0.32, 0.33"
+FFTW = "1"
+StatsBase = "0.32, 0.33, 0.34"
 julia = "1.6"
 
 [extras]

src/ToeplitzMatrices.jl

@@ -3,16 +3,28 @@ 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: parent
 import Base: ==, +, -, *, \
+import Base: AbstractMatrix
 import LinearAlgebra: Cholesky, Factorization
-import LinearAlgebra: ldiv!, factorize, lmul!, pinv, eigvals, cholesky!, cholesky, tril!, triu!, checksquare, rmul!, dot, mul!, tril, triu
+import LinearAlgebra: ldiv!, factorize, lmul!, pinv, eigvals, eigvecs, eigen, Eigen, det
+import LinearAlgebra: cholesky!, cholesky, tril!, triu!, checksquare, rmul!, dot, mul!, tril, triu
+import LinearAlgebra: istriu, istril, isdiag
 import LinearAlgebra: UpperTriangular, LowerTriangular, Symmetric, Adjoint
 import AbstractFFTs: Plan, plan_fft!
 import StatsBase
 
 export AbstractToeplitz, Toeplitz, SymmetricToeplitz, Circulant, LowerTriangularToeplitz, UpperTriangularToeplitz, TriangularToeplitz, Hankel
 export durbin, trench, levinson
 
+@static if isdefined(Base, :require_one_based_indexing)
+    const require_one_based_indexing = Base.require_one_based_indexing
+else
+    function require_one_based_indexing(A...)
+        !Base.has_offset_axes(A...) || throw(ArgumentError("offset arrays are not supported but got an array with index other than 1"))
+    end
+end
+
 include("iterativeLinearSolvers.jl")
 
 # Abstract
@@ -30,6 +42,21 @@ convert(::Type{AbstractToeplitz{T}}, A::AbstractToeplitz) where T = AbstractToep
 isconcrete(A::AbstractToeplitz) = isconcretetype(typeof(A.vc)) && isconcretetype(typeof(A.vr))
 iszero(A::AbstractToeplitz) = iszero(A.vc) && iszero(A.vr)
 
+function istril(A::AbstractToeplitz, k::Integer=0)
+    vr, vc = _vr(A), _vc(A)
+    all(iszero, @view vr[max(1, k+2):end]) && all(iszero, @view vc[2:min(-k,end)])
+end
+
+function istriu(A::AbstractToeplitz, k::Integer=0)
+    vr, vc = _vr(A), _vc(A)
+    all(iszero, @view vc[max(1, -k+2):end]) && all(iszero, @view vr[2:min(k,end)])
+end
+
+function isdiag(A::AbstractToeplitz)
+    vr, vc = _vr(A), _vc(A)
+    all(iszero, @view vr[2:end]) && all(iszero, @view vc[2:end])
+end
+
 """
     ToeplitzFactorization
 

src/hankel.jl

@@ -5,6 +5,7 @@ struct Hankel{T, V<:AbstractVector{T}, S<:DimsInteger} <: AbstractMatrix{T}
 
     function Hankel{T,V,S}(v::V, (m,n)::DimsInteger) where {T, V<:AbstractVector{T}, S<:DimsInteger}
         (m < 0 || n < 0) && throw(ArgumentError("negative size: $s"))
+        require_one_based_indexing(v)
         length(v) ≥ m+n-1 || throw(ArgumentError("inconsistency between size and number of anti-diagonals"))
         new{T,V,S}(v, (m,n))
     end
@@ -68,8 +69,7 @@ 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
+AbstractMatrix{T}(A::Hankel) where {T} = Hankel{T}(AbstractVector{T}(A.v), A.size)
 for fun in (:zero, :conj, :copy, :-, :similar, :real, :imag)
     @eval $fun(A::Hankel) = Hankel($fun(A.v), size(A))
 end

src/linearalgebra.jl

@@ -122,7 +122,7 @@ end
 function factorize(A::Toeplitz)
     T = eltype(A)
     m, n = size(A)
-    S = promote_type(T, Complex{Float32})
+    S = promote_type(float(T), Complex{Float32})
     tmp = Vector{S}(undef, m + n - 1)
     copyto!(tmp, A.vc)
     copyto!(tmp, m + 1, Iterators.reverse(A.vr), 1, n - 1)
@@ -142,7 +142,7 @@ StatsBase.levinson(A::AbstractToeplitz, B::AbstractVecOrMat) = StatsBase.levinso
 function factorize(A::SymmetricToeplitz{T}) where {T<:Number}
     vc = A.vc
     m = length(vc)
-    S = promote_type(T, Complex{Float32})
+    S = promote_type(float(T), Complex{Float32})
     tmp = Vector{S}(undef, 2 * m)
     copyto!(tmp, vc)
     @inbounds tmp[m + 1] = zero(T)
@@ -199,7 +199,7 @@ const CirculantFactorization{T, V<:AbstractVector{T}} = ToeplitzFactorization{T,
 function factorize(C::Circulant)
     T = eltype(C)
     vc = C.vc
-    S = promote_type(T, Complex{Float32})
+    S = promote_type(float(T), Complex{Float32})
     tmp = Vector{S}(undef, length(vc))
     copyto!(tmp, vc)
     dft = plan_fft!(tmp)
@@ -247,14 +247,12 @@ function ldiv!(C::CirculantFactorization, b::AbstractVector)
         ))
     end
     dft = C.dft
-    @inbounds begin
-        copyto!(tmp, b)
-        dft * tmp
-        tmp ./= vcvr_dft
-        dft \ tmp
-        T = eltype(C)
-        b .= maybereal.(T, tmp)
-    end
+    copyto!(tmp, b)
+    dft * tmp
+    tmp ./= vcvr_dft
+    dft \ tmp
+    T = eltype(C)
+    b .= maybereal.(T, tmp)
     return b
 end
 
@@ -302,6 +300,28 @@ end
 
 eigvals(C::Circulant) = eigvals(factorize(C))
 eigvals(C::CirculantFactorization) = copy(C.vcvr_dft)
+_det(C) = prod(eigvals(C))
+det(C::Circulant) = _det(C)
+det(C::Circulant{<:Real}) = real(_det(C))
+@static if VERSION <= v"1.6"
+    _cispi(x) = cis(pi*x)
+else
+    _cispi(x) = cispi(x)
+end
+function eigvecs(C::Circulant)
+    n = size(C,1)
+    M = Array{complex(float(eltype(C)))}(undef, size(C))
+    x = 2/n
+    invnorm = 1/√n
+    for CI in CartesianIndices(M)
+        k, j = Tuple(CI)
+        M[CI] = _cispi((k-1) * (j-1) * x) * invnorm
+    end
+    return M
+end
+function eigen(C::Circulant)
+    Eigen(eigvals(C), eigvecs(C))
+end
 
 sqrt(C::Circulant) = sqrt(factorize(C))
 function sqrt(C::CirculantFactorization)
@@ -394,7 +414,7 @@ function factorize(A::LowerTriangularToeplitz)
     T = eltype(A)
     v = A.v
     n = length(v)
-    S = promote_type(T, Complex{Float32})
+    S = promote_type(float(T), Complex{Float32})
     tmp = zeros(S, 2 * n - 1)
     copyto!(tmp, v)
     dft = plan_fft!(tmp)
@@ -404,7 +424,7 @@ function factorize(A::UpperTriangularToeplitz)
     T = eltype(A)
     v = A.v
     n = length(v)
-    S = promote_type(T, Complex{Float32})
+    S = promote_type(float(T), Complex{Float32})
     tmp = zeros(S, 2 * n - 1)
     tmp[1] = v[1]
     copyto!(tmp, n + 1, Iterators.reverse(v), 1, n - 1)

src/special.jl

@@ -1,50 +1,71 @@
 # special Toeplitz types that can be represented by a single vector
 # Symmetric, Circulant, LowerTriangular, UpperTriangular
+
+abstract type AbstractToeplitzSingleVector{T} <: AbstractToeplitz{T} end
+
+parent(A::AbstractToeplitzSingleVector) = A.v
+basetype(x) = basetype(typeof(x))
+
+function size(A::AbstractToeplitzSingleVector)
+    n = length(parent(A))
+    (n,n)
+end
+
+adjoint(A::AbstractToeplitzSingleVector) = transpose(conj(A))
+function zero!(A::AbstractToeplitzSingleVector)
+    fill!(parent(A), zero(eltype(A)))
+    return A
+end
+
+function lmul!(x::Number, A::AbstractToeplitzSingleVector)
+    lmul!(x,parent(A))
+    A
+end
+function rmul!(A::AbstractToeplitzSingleVector, x::Number)
+    rmul!(parent(A),x)
+    A
+end
+
+for fun in (:iszero,)
+    @eval $fun(A::AbstractToeplitzSingleVector) = $fun(parent(A))
+end
+
+AbstractToeplitz{T}(A::AbstractToeplitzSingleVector) where T = basetype(A){T}(A)
+
+(*)(scalar::Number, C::AbstractToeplitzSingleVector) = basetype(C)(scalar * parent(C))
+(*)(C::AbstractToeplitzSingleVector, scalar::Number) = basetype(C)(parent(C) * scalar)
+
+AbstractMatrix{T}(A::AbstractToeplitzSingleVector) where {T} = basetype(A){T}(AbstractVector{T}(A.v))
+
+for fun in (:zero, :conj, :copy, :-, :real, :imag)
+    @eval $fun(A::AbstractToeplitzSingleVector) = basetype(A)($fun(parent(A)))
+end
+
 for TYPE in (:SymmetricToeplitz, :Circulant, :LowerTriangularToeplitz, :UpperTriangularToeplitz)
     @eval begin
-        struct $TYPE{T, V<:AbstractVector{T}} <: AbstractToeplitz{T}
+        struct $TYPE{T, V<:AbstractVector{T}} <: AbstractToeplitzSingleVector{T}
             v::V
-        end
-        $TYPE{T}(v::V) where {T,V<:AbstractVector{T}} = $TYPE{T,V}(v)
-        $TYPE{T}(v::AbstractVector) where T = $TYPE{T}(convert(AbstractVector{T},v))
-
-        AbstractToeplitz{T}(A::$TYPE) where T = $TYPE{T}(A)
-        $TYPE{T}(A::$TYPE) where T = $TYPE{T}(convert(AbstractVector{T},A.v))
-        convert(::Type{$TYPE{T}}, A::$TYPE) where {T} = $TYPE{T}(A)
-
-        size(A::$TYPE) = (length(A.v),length(A.v))
-
-        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)
-                fill!(A.v,zero(eltype(A)))
-            else
-                A.v=zero(A.v)
+            function $TYPE{T,V}(v::V) where {T,V<:AbstractVector{T}}
+                require_one_based_indexing(v)
+                new{T,V}(v)
             end
         end
+        function $TYPE{T}(v::AbstractVector) where T
+            vT = convert(AbstractVector{T},v)
+            $TYPE{T, typeof(vT)}(vT)
+        end
+        $TYPE(v::V) where {T,V<:AbstractVector{T}} = $TYPE{T,V}(v)
+
+        basetype(::Type{T}) where {T<:$TYPE} = $TYPE
+
+        (==)(A::$TYPE, B::$TYPE) = A.v == B.v
+
+        convert(::Type{$TYPE{T}}, A::$TYPE) where {T} = A isa $TYPE{T} ? A : $TYPE{T}(A)::$TYPE{T}
 
         function copyto!(A::$TYPE,B::$TYPE)
             copyto!(A.v,B.v)
             A
         end
-        similar(A::$TYPE, T::Type) = $TYPE{T}(similar(A.v, T))
-        function lmul!(x::Number, A::$TYPE)
-            lmul!(x,A.v)
-            A
-        end
-        function rmul!(A::$TYPE, x::Number)
-            rmul!(A.v,x)
-            A
-        end
-    end
-    for fun in (:zero, :conj, :copy, :-, :real, :imag)
-        @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))
@@ -95,8 +116,12 @@ function getproperty(A::UpperTriangularToeplitz, s::Symbol)
     end
 end
 
-_circulate(v::AbstractVector) = vcat(v[1],v[end:-1:2])
-_firstnonzero(v::AbstractVector) = vcat(v[1],zero(view(v,2:lastindex(v))))
+_circulate(v::AbstractVector) = reverse(v, 2)
+function _firstnonzero(v::AbstractVector)
+    w = zero(v)
+    w[1] = v[1]
+    w
+end
 
 # transpose
 transpose(A::SymmetricToeplitz) = A
@@ -226,3 +251,6 @@ function _trisame!(A::TriangularToeplitz, k::Integer)
 end
 tril!(A::LowerTriangularToeplitz, k::Integer) = _trisame!(A,k)
 triu!(A::UpperTriangularToeplitz, k::Integer) = _trisame!(A,-k)
+
+isdiag(A::Union{Circulant, LowerTriangularToeplitz, SymmetricToeplitz}) = all(iszero, @view _vc(A)[2:end])
+isdiag(A::UpperTriangularToeplitz) = all(iszero, @view _vr(A)[2:end])

src/toeplitz.jl

@@ -9,6 +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}}
+        require_one_based_indexing(vr, vc)
         if first(vc) != first(vr)
             error("First element of the vectors must be the same")
         end
@@ -104,10 +105,9 @@ 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]
+function AbstractMatrix{T}(A::AbstractToeplitz) where {T}
+    vc = AbstractVector{T}(_vc(A))
+    vr = AbstractVector{T}(_vr(A))
     Toeplitz{T}(vc,vr)
 end
 for fun in (:zero, :conj, :copy, :-, :real, :imag)