From fde025f2fe3643e8673a098fed02bc0d804dc7ed Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Fri, 4 Aug 2023 13:28:21 +0400 Subject: [PATCH 01/25] Split y coordinate computation in `paduapoints` (#226) * split y coordinate computation in paduapoints * Add inbounds and tests * Bump version to v0.15.7 --- Project.toml | 2 +- src/PaduaTransform.jl | 41 +++++++++++++++++++++++++++++++++-------- test/paduatests.jl | 13 +++++++++++++ 3 files changed, 47 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index 34c7af00..00bb0171 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.15.6" +version = "0.15.7" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/src/PaduaTransform.jl b/src/PaduaTransform.jl index 6188e74d..c7670924 100644 --- a/src/PaduaTransform.jl +++ b/src/PaduaTransform.jl @@ -209,21 +209,46 @@ function paduapoints(::Type{T}, n::Integer) where T MM=Matrix{T}(undef,N,2) m=0 delta=0 - NN=fld(n+2,2) - @inbounds for k=n:-1:0 - if isodd(n)>0 - delta=mod(k,2) + NN=div(n,2)+1 + # x coordinates + for k=n:-1:0 + if isodd(n) + delta = Int(isodd(k)) end + x = -cospi(T(k)/n) @inbounds for j=NN+delta:-1:1 m+=1 - MM[m,1]=sinpi(T(k)/n-T(0.5)) - if isodd(n-k)>0 - MM[m,2]=sinpi((2j-one(T))/(n+1)-T(0.5)) + MM[m,1]=x + end + end + # y coordinates + # populate the first two sets, and copy the rest + m=0 + for k=n:-1:n-1 + if isodd(n) + delta = Int(isodd(k)) + end + for j=NN+delta:-1:1 + m+=1 + @inbounds if isodd(n-k) + MM[m,2]=-cospi((2j-one(T))/(n+1)) else - MM[m,2]=sinpi(T(2j-2)/(n+1)-T(0.5)) + MM[m,2]=-cospi(T(2j-2)/(n+1)) end end end + m += 1 + # number of y coordinates between k=n and k=n-2 + Ny_shift = 2NN+isodd(n) + for k in n-2:-1:0 + if isodd(n) + delta = Int(isodd(k)) + end + for j in range(m, length=NN+delta) + @inbounds MM[j,2] = MM[j-Ny_shift,2] + end + m += NN+delta + end return MM end diff --git a/test/paduatests.jl b/test/paduatests.jl index cc46d462..c82dc579 100644 --- a/test/paduatests.jl +++ b/test/paduatests.jl @@ -53,4 +53,17 @@ using FastTransforms, Test g_l=paduaeval(g_xy,x,y,l,Val{false}) @test f_xy(x,y) ≈ f_m @test g_xy(x,y) ≈ g_l + + # odd n + m=135 + l=85 + f_m=paduaeval(f_xy,x,y,m,Val{true}) + g_l=paduaeval(g_xy,x,y,l,Val{true}) + @test f_xy(x,y) ≈ f_m + @test g_xy(x,y) ≈ g_l + + f_m=paduaeval(f_xy,x,y,m,Val{false}) + g_l=paduaeval(g_xy,x,y,l,Val{false}) + @test f_xy(x,y) ≈ f_m + @test g_xy(x,y) ≈ g_l end From 6fd000cf130b624b9a65a55127dafd03fe7adf37 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Tue, 3 Oct 2023 11:34:41 +0100 Subject: [PATCH 02/25] CompatHelper: bump compat for FastGaussQuadrature to 1, (keep existing compat) (#227) * CompatHelper: bump compat for FastGaussQuadrature to 1, (keep existing compat) * v0.15.8 --------- Co-authored-by: CompatHelper Julia Co-authored-by: Sheehan Olver --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 00bb0171..581a185c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.15.7" +version = "0.15.8" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" @@ -19,7 +19,7 @@ ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24" [compat] AbstractFFTs = "1.0" FFTW = "1.7" -FastGaussQuadrature = "0.4, 0.5" +FastGaussQuadrature = "0.4, 0.5, 1" FastTransforms_jll = "0.6.2" FillArrays = "0.9, 0.10, 0.11, 0.12, 0.13, 1" GenericFFT = "0.1" From 50282780c8a377ae700109a5e76df7e81b3c14f1 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 21 Oct 2023 17:36:24 +0100 Subject: [PATCH 03/25] matrix th_ultra2ultra, th_jac2jac (#228) * matrix th_ultra2ultra * leg2chebu tests * jac2jac matrix * support composing plans in Toeplitz-Hankel * Mape type stable via vector plans * prepare for integer increments * integer increments with vectors * fix complexity of integer increment * JacJacPlan remembers dims * matrix raising * Fix Chebyshev transform * higher th_ultra2ultra * no-alloc jacobi integer step * fix tests * fix nearest integer --- Project.toml | 2 +- src/toeplitzhankel.jl | 400 ++++++++++++++++++++++++++++++++++-- test/toeplitzhankeltests.jl | 74 ++++++- 3 files changed, 455 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index 581a185c..8a132703 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.15.8" +version = "0.15.9" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/src/toeplitzhankel.jl b/src/toeplitzhankel.jl index ff9acc52..eb11994e 100644 --- a/src/toeplitzhankel.jl +++ b/src/toeplitzhankel.jl @@ -1,7 +1,20 @@ """ -Store a diagonally-scaled Toeplitz∘Hankel matrix: +Represent a scaled Toeplitz∘Hankel matrix: + DL(T∘H)DR -where the Hankel matrix `H` is non-negative definite. This allows a Cholesky decomposition in 𝒪(K²N) operations and 𝒪(KN) storage, K = log N log ɛ⁻¹. + +where the Hankel matrix `H` is non-negative definite, via + + ∑_{k=1}^r Diagonal(L[:,k])*T*Diagonal(R[:,k]) + +where `L` and `R` are determined by doing a rank-r pivoted Cholesky decomposition of `H`, which in low rank form is + + H ≈ ∑_{k=1}^r C[:,k]C[:,k]' + +so that `L[:,k] = DL*C[:,k]` and `R[:,k] = DR*C[:,k]`. + +This allows a Cholesky decomposition in 𝒪(K²N) operations and 𝒪(KN) storage, K = log N log ɛ⁻¹. +The tuple storage allows plans applied to each dimension. """ struct ToeplitzHankelPlan{S, N, M, N1, TP<:ToeplitzPlan{S,N1}} <: Plan{S} T::NTuple{M,TP} @@ -209,7 +222,7 @@ function _leg2chebuTH_TLC(::Type{S}, mn, d) where {S} t[1:2:end] = λ[1:2:n]./(((1:2:n).-2)) h = λ./((1:2n-1).+1) C = hankel_partialchol(h) - T = plan_uppertoeplitz!(-2t/π, (length(t), size(C,2)), 1) + T = plan_uppertoeplitz!(-2t/π, (mn..., size(C,2)), d) (T, (1:n) .* C, C) end @@ -229,6 +242,10 @@ for f in (:leg2cheb, :leg2chebu) end end +### +# th_cheb2leg +### + _sub_dim_by_one(d) = () _sub_dim_by_one(d, m, n...) = (isone(d) ? m-1 : m, _sub_dim_by_one(d-1, n...)...) @@ -257,7 +274,38 @@ function plan_th_cheb2leg!(::Type{S}, mn::NTuple{2,Int}, dims::NTuple{2,Int}) wh ChebyshevToLegendrePlanTH(ToeplitzHankelPlan((T1,T2), (L1,L2), (C1,C2), dims)) end -function plan_th_ultra2ultra!(::Type{S}, (n,)::Tuple{Int}, λ₁, λ₂) where {S} + +### +# th_ultra2ultra +### + +# The second case handles zero +isapproxinteger(::Integer) = true +isapproxinteger(x) = isinteger(x) || x ≈ round(Int,x) || x+1 ≈ round(Int,x+1) + + +_nearest_jacobi_par(α, γ) = isapproxinteger(α-γ) ? α : round(Int,α,RoundDown) + mod(γ,1) + + +struct Ultra2UltraPlanTH{T, Plans, Dims} <: Plan{T} + plans::Plans + λ₁::T + λ₂::T + dims::Dims +end + +function *(P::Ultra2UltraPlanTH, A::AbstractArray) + ret = A + for p in P.plans + ret = p*ret + end + c = _nearest_jacobi_par(P.λ₁, P.λ₂) + + _ultra2ultra_integerinc!(ret, c, P.λ₂, P.dims) +end + +function _ultra2ultraTH_TLC(::Type{S}, mn, λ₁, λ₂, d) where {S} + n = mn[d] @assert abs(λ₁-λ₂) < 1 S̃ = real(S) DL = (zero(S̃):n-one(S̃)) .+ λ₂ @@ -267,8 +315,151 @@ function plan_th_ultra2ultra!(::Type{S}, (n,)::Tuple{Int}, λ₁, λ₂) where { h = Λ.(jk,λ₁,λ₂+one(S̃)) lmul!(gamma(λ₂)/gamma(λ₁),h) C = hankel_partialchol(h) - T = plan_uppertoeplitz!(lmul!(inv(gamma(λ₁-λ₂)),t), (length(t), size(C,2)), 1) - ToeplitzHankelPlan(T, DL .* C, C) + T = plan_uppertoeplitz!(lmul!(inv(gamma(λ₁-λ₂)),t), (mn..., size(C,2)), d) + T, DL .* C, C +end + +_good_plan_th_ultra2ultra!(::Type{S}, mn, λ₁, λ₂, dims::Int) where S = ToeplitzHankelPlan(_ultra2ultraTH_TLC(S, mn, λ₁, λ₂, dims)..., dims) + +function _good_plan_th_ultra2ultra!(::Type{S}, mn::NTuple{2,Int}, λ₁, λ₂, dims::NTuple{2,Int}) where S + T1,L1,C1 = _ultra2ultraTH_TLC(S, mn, λ₁, λ₂, 1) + T2,L2,C2 = _ultra2ultraTH_TLC(S, mn, λ₁, λ₂, 2) + ToeplitzHankelPlan((T1,T2), (L1,L2), (C1,C2), dims) +end + + + +function plan_th_ultra2ultra!(::Type{S}, mn, λ₁, λ₂, dims) where {S} + c = _nearest_jacobi_par(λ₁, λ₂) + + if isapproxinteger(λ₂ - λ₁) + # TODO: don't make extra plan + plans = typeof(_good_plan_th_ultra2ultra!(S, mn, λ₁+0.1, λ₂, dims))[] + else + plans = [_good_plan_th_ultra2ultra!(S, mn, λ₁, c, dims)] + end + + Ultra2UltraPlanTH(plans, λ₁, λ₂, dims) +end + +function _ultra_raise!(B, λ) + m, n = size(B, 1), size(B, 2) + + if m > 1 + @inbounds for j = 1:n + for i = 1:m-2 + Bij = λ / (i+λ-1) * B[i,j] + Bij += -λ / (i+λ+1) * B[i+2,j] + B[i,j] = Bij + end + B[m-1,j] = λ / (m+λ-2)*B[m-1,j] + B[m,j] = λ / (m+λ-1)*B[m,j] + end + end + B +end + +function _ultra_lower!(B, λ) + m, n = size(B, 1), size(B, 2) + + if m > 1 + @inbounds for j = 1:n + B[m,j] = (m+λ-1)/λ * B[m,j] + B[m-1,j] = (m+λ-2)/λ *B[m-1,j] + for i = m-2:-1:1 + Bij = B[i,j] + λ / (i+λ+1) * B[i+2,j] + B[i,j] = (i+λ-1)/λ * Bij + end + end + end + B +end + + + +function _ultra_raise!(x, λ, dims) + for d in dims + if d == 1 + _ultra_raise!(x, λ) + else + _ultra_raise!(x', λ) + end + end + x +end + +function _ultra_lower!(x, λ, dims) + for d in dims + if d == 1 + _ultra_lower!(x, λ-1) + else + _ultra_lower!(x', λ-1) + end + end + x +end + +function _ultra2ultra_integerinc!(x, λ₁, λ₂, dims) + while !(λ₁ ≈ λ₂) + if λ₂ > λ₁ + _ultra_raise!(x, λ₁, dims) + λ₁ += 1 + else + _ultra_lower!(x, λ₁, dims) + λ₁ -= 1 + end + end + x +end + +### +# th_jac2jac +### + + +function _lmul!(A::Bidiagonal, B::AbstractVecOrMat) + @assert A.uplo == 'U' + + m, n = size(B, 1), size(B, 2) + if m != size(A, 1) + throw(DimensionMismatch("right hand side B needs first dimension of size $(size(A,1)), has size $m")) + end + @inbounds for j = 1:n + for i = 1:m-1 + Bij = A.dv[i]*B[i,j] + Bij += A.ev[i]*B[i+1,j] + B[i,j] = Bij + end + B[m,j] = A.dv[m]*B[m,j] + end + B +end + +struct Jac2JacPlanTH{T, Plans, Dims} <: Plan{T} + plans::Plans + α::T + β::T + γ::T + δ::T + dims::Dims +end + +Jac2JacPlanTH(plans, α, β, γ, δ, dims) = Jac2JacPlanTH(plans, promote(α, β, γ, δ)..., dims) + +function *(P::Jac2JacPlanTH, A::AbstractArray) + if P.α + P.β ≤ -1 + _jacobi_raise_a!(A, P.α, P.β) + c,d = _nearest_jacobi_par(P.α+1, P.γ), _nearest_jacobi_par(P.β, P.δ) + else + c,d = _nearest_jacobi_par(P.α, P.γ), _nearest_jacobi_par(P.β, P.δ) + end + + ret = A + for p in P.plans + ret = p*ret + end + + _jac2jac_integerinc!(ret, c, d, P.γ, P.δ, P.dims) end function alternatesign!(v) @@ -278,32 +469,206 @@ function alternatesign!(v) v end -function plan_th_jac2jac!(::Type{S}, (n,), α, β, γ, δ) where {S} +function _jac2jacTH_TLC(::Type{S}, mn, α, β, γ, δ, d) where {S} + n = mn[d] + @assert α+β > -1 if β == δ @assert abs(α-γ) < 1 - @assert α+β > -1 jk = 0:n-1 DL = (2jk .+ γ .+ β .+ 1).*Λ.(jk,γ+β+1,β+1) t = convert(AbstractVector{S}, Λ.(jk, α-γ,1)) h = Λ.(0:2n-2,α+β+1,γ+β+2) DR = Λ.(jk,β+1,α+β+1)./gamma(α-γ) C = hankel_partialchol(h) - T = plan_uppertoeplitz!(t, (length(t), size(C,2)), 1) + T = plan_uppertoeplitz!(t, (mn..., size(C,2)), d) elseif α == γ + @assert abs(β-δ) < 1 jk = 0:n-1 DL = (2jk .+ δ .+ α .+ 1).*Λ.(jk,δ+α+1,α+1) h = Λ.(0:2n-2,α+β+1,δ+α+2) DR = Λ.(jk,α+1,α+β+1)./gamma(β-δ) t = alternatesign!(convert(AbstractVector{S}, Λ.(jk,β-δ,1))) C = hankel_partialchol(h) - T = plan_uppertoeplitz!(t, (length(t), size(C,2)), 1) + T = plan_uppertoeplitz!(t, (mn..., size(C,2)), d) else throw(ArgumentError("Cannot create Toeplitz dot Hankel, use a sequence of plans.")) end - ToeplitzHankelPlan(T, DL .* C, DR .* C) + (T, DL .* C, DR .* C) +end + +_good_plan_th_jac2jac!(::Type{S}, mn, α, β, γ, δ, dims::Int) where S = ToeplitzHankelPlan(_jac2jacTH_TLC(S, mn, α, β, γ, δ, dims)..., dims) + +function _good_plan_th_jac2jac!(::Type{S}, mn::NTuple{2,Int}, α, β, γ, δ, dims::NTuple{2,Int}) where S + T1,L1,C1 = _jac2jacTH_TLC(S, mn, α, β, γ, δ, 1) + T2,L2,C2 = _jac2jacTH_TLC(S, mn, α, β, γ, δ, 2) + ToeplitzHankelPlan((T1,T2), (L1,L2), (C1,C2), dims) +end + + + +function plan_th_jac2jac!(::Type{S}, mn, α, β, γ, δ, dims) where {S} + if α + β ≤ -1 + c,d = _nearest_jacobi_par(α+1, γ), _nearest_jacobi_par(β, δ) + else + c,d = _nearest_jacobi_par(α, γ), _nearest_jacobi_par(β, δ) + end + + if isapproxinteger(β - δ) && isapproxinteger(α-γ) + # TODO: don't make extra plan + plans = typeof(_good_plan_th_jac2jac!(S, mn, α+0.1, β, α, β, dims))[] + elseif isapproxinteger(α - γ) || isapproxinteger(β - δ) + if α + β ≤ -1 + # avoid degenerecies + plans = [_good_plan_th_jac2jac!(S, mn, α+1, β, c, d, dims)] + else + plans = [_good_plan_th_jac2jac!(S, mn, α, β, c, d, dims)] + end + else + if α + β ≤ -1 + plans = [_good_plan_th_jac2jac!(S, mn, α+1, β, α+1, d, dims), _good_plan_th_jac2jac!(S, mn, α+1, d, c, d, dims)] + else + plans = [_good_plan_th_jac2jac!(S, mn, α, β, α, d, dims), _good_plan_th_jac2jac!(S, mn, α, d, c, d, dims)] + end + end + + Jac2JacPlanTH(plans, α, β, γ, δ, dims) +end + + +function _jacobi_raise_a!(B, a, b) + m, n = size(B, 1), size(B, 2) + if m > 1 + @inbounds for j = 1:n + B[1,j] = B[1,j] - (1+b) / (a+b+3) * B[2,j] + for i = 2:m-1 + B[i,j] = (i+a+b)/(a+b-1+2i) * B[i,j] - (i+b) / (a+b+2i+1) * B[i+1,j] + end + B[m,j] = (m+a+b)/(a+b-1+2m)*B[m,j] + end + end + B +end + +function _jacobi_lower_a!(B, a, b) + m, n = size(B, 1), size(B, 2) + + if m > 1 + @inbounds for j = 1:n + B[m,j] = (a+b-1+2m)/(m+a+b) * B[m,j] + for i = m-1:-1:2 + Bij = B[i,j] + (i+b) / (a+b+2i+1) * B[i+1,j] + B[i,j] = (a+b-1+2i)/(i+a+b) * Bij + end + B[1,j] = B[1,j] + (1+b) / (a+b+3) * B[2,j] + end + end + B +end + + + +function _jacobi_raise_b!(B, a, b) + m, n = size(B, 1), size(B, 2) + if m > 1 + @inbounds for j = 1:n + B[1,j] = B[1,j] + (1+a) / (a+b+3) * B[2,j] + + for i = 2:m-1 + B[i,j] = (i+a+b)/(a+b-1+2i) * B[i,j] + (i+a) / (a+b+2i+1) * B[i+1,j] + end + B[m,j] = (m+a+b)/(a+b-1+2m)*B[m,j] + end + end + B +end + +function _jacobi_lower_b!(B, a, b) + m, n = size(B, 1), size(B, 2) + + if m > 1 + @inbounds for j = 1:n + B[m,j] = (a+b-1+2m)/(m+a+b) * B[m,j] + for i = m-1:-1:2 + Bij = B[i,j] - (i+a) / (a+b+2i+1) * B[i+1,j] + B[i,j] = (a+b-1+2i)/(i+a+b) * Bij + end + B[1,j] = B[1,j] - (1+a) / (a+b+3) * B[2,j] + end + end + B +end + + + +function _jacobi_raise_b!(x, α, β, dims) + for d in dims + if d == 1 + _jacobi_raise_b!(x, α, β) + else + _jacobi_raise_b!(x', α, β) + end + end + x +end +function _jacobi_raise_a!(x, α, β, dims) + for d in dims + if d == 1 + _jacobi_raise_a!(x, α, β) + else + _jacobi_raise_a!(x', α, β) + end + end + x +end + +function _jacobi_lower_b!(x, α, β, dims) + for d in dims + if d == 1 + _jacobi_lower_b!(x, α, β-1) + else + _jacobi_lower_b!(x', α, β-1) + end + end + x +end +function _jacobi_lower_a!(x, α, β, dims) + for d in dims + if d == 1 + _jacobi_lower_a!(x, α-1, β) + else + _jacobi_lower_a!(x', α-1, β) + end + end + x +end + + +function _jac2jac_integerinc!(x, α, β, γ, δ, dims) + while !(α ≈ γ && β ≈ δ) + if !(δ ≈ β) && δ > β + _jacobi_raise_b!(x, α, β, dims) + β += 1 + elseif !(δ ≈ β) && δ < β + _jacobi_lower_b!(x, α, β, dims) + β -= 1 + elseif !(γ ≈ α) && γ > α + _jacobi_raise_a!(x, α, β, dims) + α += 1 + else + @assert γ < α + _jacobi_lower_a!(x, α, β, dims) + α -= 1 + end + end + x end + +### +# other routines +### + for f in (:th_leg2cheb, :th_cheb2leg, :th_leg2chebu) plan = Symbol("plan_", f, "!") @eval begin @@ -315,5 +680,14 @@ for f in (:th_leg2cheb, :th_cheb2leg, :th_leg2chebu) end end -th_ultra2ultra(v, λ₁, λ₂, dims...) = plan_th_ultra2ultra!(eltype(v),size(v),λ₁,λ₂, dims...)*copy(v) -th_jac2jac(v, α, β, γ, δ, dims...) = plan_th_jac2jac!(eltype(v),size(v),α,β,γ,δ, dims...)*copy(v) \ No newline at end of file +plan_th_ultra2ultra!(::Type{S}, mn::NTuple{N,Int}, λ₁, λ₂, dims::UnitRange) where {N,S} = plan_th_ultra2ultra!(S, mn, λ₁, λ₂, tuple(dims...)) +plan_th_ultra2ultra!(::Type{S}, mn::Tuple{Int}, λ₁, λ₂, dims::Tuple{Int}=(1,)) where {S} = plan_th_ultra2ultra!(S, mn, λ₁, λ₂, dims...) +plan_th_ultra2ultra!(::Type{S}, (m,n)::NTuple{2,Int}, λ₁, λ₂) where {S} = plan_th_ultra2ultra!(S, (m,n), λ₁, λ₂, (1,2)) +plan_th_ultra2ultra!(arr::AbstractArray{T}, λ₁, λ₂, dims...) where T = plan_th_ultra2ultra!(T, size(arr), λ₁, λ₂, dims...) +th_ultra2ultra(v, λ₁, λ₂, dims...) = plan_th_ultra2ultra!(eltype(v), size(v), λ₁, λ₂, dims...)*copy(v) + +plan_th_jac2jac!(::Type{S}, mn::NTuple{N,Int}, α, β, γ, δ, dims::UnitRange) where {N,S} = plan_th_jac2jac!(S, mn, α, β, γ, δ, tuple(dims...)) +plan_th_jac2jac!(::Type{S}, mn::Tuple{Int}, α, β, γ, δ, dims::Tuple{Int}=(1,)) where {S} = plan_th_jac2jac!(S, mn, α, β, γ, δ, dims...) +plan_th_jac2jac!(::Type{S}, (m,n)::NTuple{2,Int}, α, β, γ, δ) where {S} = plan_th_jac2jac!(S, (m,n), α, β, γ, δ, (1,2)) +plan_th_jac2jac!(arr::AbstractArray{T}, α, β, γ, δ, dims...) where T = plan_th_jac2jac!(T, size(arr), α, β, γ, δ, dims...) +th_jac2jac(v, α, β, γ, δ, dims...) = plan_th_jac2jac!(eltype(v), size(v), α, β, γ, δ, dims...)*copy(v) \ No newline at end of file diff --git a/test/toeplitzhankeltests.jl b/test/toeplitzhankeltests.jl index a3f8d36c..ce1d2f69 100644 --- a/test/toeplitzhankeltests.jl +++ b/test/toeplitzhankeltests.jl @@ -1,7 +1,7 @@ using FastTransforms, Test -import FastTransforms: th_leg2cheb, th_cheb2leg, th_ultra2ultra,th_jac2jac, th_leg2chebu, +import FastTransforms: th_leg2cheb, th_cheb2leg, th_leg2chebu, th_ultra2ultra,th_jac2jac, th_leg2chebu, lib_leg2cheb, lib_cheb2leg, lib_ultra2ultra, lib_jac2jac, - plan_th_cheb2leg!, plan_th_leg2cheb! + plan_th_cheb2leg!, plan_th_leg2chebu!, plan_th_leg2cheb!, plan_th_ultra2ultra!, plan_th_jac2jac! @testset "ToeplitzHankel" begin for x in ([1.0], [1.0,2,3,4,5], [1.0+im,2-3im,3+4im,4-5im,5+10im], collect(1.0:1000)) @@ -9,17 +9,38 @@ import FastTransforms: th_leg2cheb, th_cheb2leg, th_ultra2ultra,th_jac2jac, th_l @test th_cheb2leg(x) ≈ lib_cheb2leg(x) @test th_leg2chebu(x) ≈ lib_ultra2ultra(x, 0.5, 1.0) @test th_ultra2ultra(x,0.1, 0.2) ≈ lib_ultra2ultra(x, 0.1, 0.2) - @test th_jac2jac(x,0.1, 0.2,0.1,0.4) ≈ lib_jac2jac(x, 0.1, 0.2,0.1,0.4) + @test th_ultra2ultra(x,1, 2) ≈ lib_ultra2ultra(x, 1, 2) + @test th_ultra2ultra(x,0.1, 2.2) ≈ lib_ultra2ultra(x, 0.1, 2.2) + @test th_ultra2ultra(x, 2.2, 0.1) ≈ lib_ultra2ultra(x, 2.2, 0.1) + @test @inferred(th_jac2jac(x,0.1, 0.2,0.1,0.4)) ≈ lib_jac2jac(x, 0.1, 0.2,0.1,0.4) @test th_jac2jac(x,0.1, 0.2,0.3,0.2) ≈ lib_jac2jac(x, 0.1, 0.2,0.3,0.2) + @test th_jac2jac(x,0.1, 0.2,0.3,0.4) ≈ lib_jac2jac(x, 0.1, 0.2,0.3,0.4) + @test @inferred(th_jac2jac(x,0.1, 0.2,1.3,0.4)) ≈ lib_jac2jac(x, 0.1, 0.2,1.3,0.4) + @test th_jac2jac(x,0.1, 0.2,1.3,2.4) ≈ lib_jac2jac(x, 0.1, 0.2,1.3,2.4) + @test th_jac2jac(x,1.3,2.4, 0.1, 0.2) ≈ lib_jac2jac(x,1.3,2.4, 0.1, 0.2) + @test th_jac2jac(x,1.3, 1.2,-0.1,-0.2) ≈ lib_jac2jac(x, 1.3, 1.2,-0.1,-0.2) + @test @inferred(th_jac2jac(x,-0.5, -0.5, -0.5,-0.5)) ≈ lib_jac2jac(x, -0.5, -0.5, -0.5,-0.5) + @test th_jac2jac(x,-0.5, -0.5, 0.5,0.5) ≈ lib_jac2jac(x, -0.5, -0.5, 0.5,0.5) + @test th_jac2jac(x,0.5,0.5,-0.5, -0.5) ≈ lib_jac2jac(x, 0.5,0.5,-0.5, -0.5) + @test th_jac2jac(x,-0.5, -0.5, 0.5,-0.5) ≈ lib_jac2jac(x, -0.5, -0.5, 0.5,-0.5) + @test th_jac2jac(x, -1/2,-1/2,1/2,0) ≈ lib_jac2jac(x, -1/2,-1/2,1/2,0) + @test th_jac2jac(x, -1/2,-1/2,0,1/2) ≈ lib_jac2jac(x, -1/2,-1/2,0,1/2) + @test th_jac2jac(x, -3/4,-3/4,0,3/4) ≈ lib_jac2jac(x, -3/4,-3/4,0,3/4) + @test th_jac2jac(x,0, 0, 5, 5) ≈ lib_jac2jac(x, 0, 0, 5, 5) + if length(x) < 10 + @test th_jac2jac(x, 5, 5, 0, 0) ≈ lib_jac2jac(x, 5, 5, 0, 0) + end - - @test th_cheb2leg(th_leg2cheb(x)) ≈ x atol=1E-9 - @test th_leg2cheb(th_cheb2leg(x)) ≈ x atol=1E-10 + @test th_cheb2leg(th_leg2cheb(x)) ≈ x + @test th_leg2cheb(th_cheb2leg(x)) ≈ x + @test th_ultra2ultra(th_ultra2ultra(x, 0.1, 0.6), 0.6, 0.1) ≈ x + @test th_jac2jac(th_jac2jac(x, 0.1, 0.6, 0.1, 0.8), 0.1, 0.8, 0.1, 0.6) ≈ x + @test th_jac2jac(th_jac2jac(x, 0.1, 0.6, 0.2, 0.8), 0.2, 0.8, 0.1, 0.6) ≈ x end for X in (randn(5,4), randn(5,4) + im*randn(5,4)) @test th_leg2cheb(X, 1) ≈ hcat([leg2cheb(X[:,j]) for j=1:size(X,2)]...) - @test_broken th_leg2cheb(X, 1) ≈ leg2cheb(X, 1) + @test_broken th_leg2cheb(X, 1) ≈ leg2cheb(X, 1) # matrices not supported in FastTransforms @test th_leg2cheb(X, 2) ≈ vcat([permutedims(leg2cheb(X[k,:])) for k=1:size(X,1)]...) @test_broken th_leg2cheb(X, 2) ≈ leg2cheb(X, 2) @test th_leg2cheb(X) ≈ th_leg2cheb(th_leg2cheb(X, 1), 2) @@ -33,6 +54,45 @@ import FastTransforms: th_leg2cheb, th_cheb2leg, th_ultra2ultra,th_jac2jac, th_l @test th_leg2cheb(X) == plan_th_leg2cheb!(X, 1:2)*copy(X) @test th_leg2cheb(th_cheb2leg(X)) ≈ X + + @test th_leg2chebu(X, 1) ≈ hcat([ultra2ultra(X[:,j], 0.5, 1.0) for j=1:size(X,2)]...) + @test th_leg2chebu(X, 2) ≈ vcat([permutedims(ultra2ultra(X[k,:], 0.5, 1.0)) for k=1:size(X,1)]...) + @test th_leg2chebu(X) ≈ th_leg2chebu(th_leg2chebu(X, 1), 2) + + @test th_leg2chebu(X) == plan_th_leg2chebu!(X, 1:2)*copy(X) + + @test th_ultra2ultra(X, 0.1, 0.6, 1) ≈ hcat([ultra2ultra(X[:,j], 0.1, 0.6) for j=1:size(X,2)]...) + @test th_ultra2ultra(X, 0.1, 0.6, 2) ≈ vcat([permutedims(ultra2ultra(X[k,:], 0.1, 0.6)) for k=1:size(X,1)]...) + @test th_ultra2ultra(X, 0.1, 0.6) ≈ th_ultra2ultra(th_ultra2ultra(X, 0.1, 0.6, 1), 0.1, 0.6, 2) + + @test th_ultra2ultra(X, 0.1, 2.6, 1) ≈ hcat([ultra2ultra(X[:,j], 0.1, 2.6) for j=1:size(X,2)]...) + @test th_ultra2ultra(X, 0.1, 2.6, 2) ≈ vcat([permutedims(ultra2ultra(X[k,:], 0.1, 2.6)) for k=1:size(X,1)]...) + @test th_ultra2ultra(X, 0.1, 2.6) ≈ th_ultra2ultra(th_ultra2ultra(X, 0.1, 2.6, 1), 0.1, 2.6, 2) + + @test th_ultra2ultra(X, 2.6, 0.1, 1) ≈ hcat([ultra2ultra(X[:,j], 2.6, 0.1) for j=1:size(X,2)]...) + @test th_ultra2ultra(X, 2.6, 0.1, 2) ≈ vcat([permutedims(ultra2ultra(X[k,:], 2.6, 0.1)) for k=1:size(X,1)]...) + @test th_ultra2ultra(X, 2.6, 0.1) ≈ th_ultra2ultra(th_ultra2ultra(X, 2.6, 0.1, 1), 2.6, 0.1, 2) + + @test th_ultra2ultra(X, 0.1, 0.6) == plan_th_ultra2ultra!(X, 0.1, 0.6, 1:2)*copy(X) + @test th_ultra2ultra(X, 0.1, 0.6) == plan_th_ultra2ultra!(X, 0.1, 0.6, 1:2)*copy(X) + + @test th_ultra2ultra(th_ultra2ultra(X, 0.1, 0.6), 0.6, 0.1) ≈ X + + @test th_jac2jac(X, 0.1, 0.6, 0.1, 0.8, 1) ≈ hcat([jac2jac(X[:,j], 0.1, 0.6, 0.1, 0.8) for j=1:size(X,2)]...) + @test th_jac2jac(X, 0.1, 0.6, 0.1, 0.8, 2) ≈ vcat([permutedims(jac2jac(X[k,:], 0.1, 0.6, 0.1, 0.8)) for k=1:size(X,1)]...) + @test th_jac2jac(X, 0.1, 0.6, 0.1, 0.8) ≈ th_jac2jac(th_jac2jac(X, 0.1, 0.6, 0.1, 0.8, 1), 0.1, 0.6, 0.1, 0.8, 2) + + @test th_jac2jac(X, 0.1, 0.6, 0.2, 0.8, 1) ≈ hcat([jac2jac(X[:,j], 0.1, 0.6, 0.2, 0.8) for j=1:size(X,2)]...) + @test th_jac2jac(X, 0.1, 0.6, 0.2, 0.8, 2) ≈ vcat([permutedims(jac2jac(X[k,:], 0.1, 0.6, 0.2, 0.8)) for k=1:size(X,1)]...) + + @test th_jac2jac(X, 0.1, 0.6, 0.1, 0.8) == plan_th_jac2jac!(X, 0.1, 0.6, 0.1, 0.8, 1:2)*copy(X) + @test th_jac2jac(X, 0.1, 0.6, 0.1, 0.8) == plan_th_jac2jac!(X, 0.1, 0.6, 0.1, 0.8, 1:2)*copy(X) + + @test th_jac2jac(th_jac2jac(X, 0.1, 0.6, 0.1, 0.8), 0.1, 0.8, 0.1, 0.6) ≈ X + + @test th_jac2jac(X, 0.1, 0.6, 3.1, 2.8, 1) ≈ hcat([jac2jac(X[:,j], 0.1, 0.6, 3.1, 2.8) for j=1:size(X,2)]...) + @test th_jac2jac(X, 0.1, 0.6, 3.1, 2.8, 2) ≈ vcat([permutedims(jac2jac(X[k,:], 0.1, 0.6, 3.1, 2.8)) for k=1:size(X,1)]...) + @test th_jac2jac(X, 0.1, 0.6, 3.1, 2.8) ≈ th_jac2jac(th_jac2jac(X, 0.1, 0.6, 3.1, 2.8, 1), 0.1, 0.6, 3.1, 2.8, 2) end @testset "BigFloat" begin From 2dcabebd3b31a59b3ece09e390e127a279a941c8 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 21 Oct 2023 21:41:30 +0100 Subject: [PATCH 04/25] Add th_cheb2jac and th_jac2cheb (#229) * Add Cheb2Jac and Jac2Cheb * seed random nums * Update toeplitzhankeltests.jl --- Project.toml | 9 ++++- src/toeplitzhankel.jl | 72 +++++++++++++++++++++++++++++++++++-- test/toeplitzhankeltests.jl | 24 +++++++++++-- 3 files changed, 99 insertions(+), 6 deletions(-) diff --git a/Project.toml b/Project.toml index 8a132703..3a97bbd3 100644 --- a/Project.toml +++ b/Project.toml @@ -13,7 +13,6 @@ Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24" [compat] @@ -27,3 +26,11 @@ Reexport = "0.2, 1.0" SpecialFunctions = "0.10, 1, 2" ToeplitzMatrices = "0.7.1, 0.8" julia = "1.7" + + +[extras] +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test", "Random"] diff --git a/src/toeplitzhankel.jl b/src/toeplitzhankel.jl index eb11994e..3bd09d55 100644 --- a/src/toeplitzhankel.jl +++ b/src/toeplitzhankel.jl @@ -448,7 +448,7 @@ Jac2JacPlanTH(plans, α, β, γ, δ, dims) = Jac2JacPlanTH(plans, promote(α, β function *(P::Jac2JacPlanTH, A::AbstractArray) if P.α + P.β ≤ -1 - _jacobi_raise_a!(A, P.α, P.β) + _jacobi_raise_a!(A, P.α, P.β, P.dims) c,d = _nearest_jacobi_par(P.α+1, P.γ), _nearest_jacobi_par(P.β, P.δ) else c,d = _nearest_jacobi_par(P.α, P.γ), _nearest_jacobi_par(P.β, P.δ) @@ -690,4 +690,72 @@ plan_th_jac2jac!(::Type{S}, mn::NTuple{N,Int}, α, β, γ, δ, dims::UnitRange) plan_th_jac2jac!(::Type{S}, mn::Tuple{Int}, α, β, γ, δ, dims::Tuple{Int}=(1,)) where {S} = plan_th_jac2jac!(S, mn, α, β, γ, δ, dims...) plan_th_jac2jac!(::Type{S}, (m,n)::NTuple{2,Int}, α, β, γ, δ) where {S} = plan_th_jac2jac!(S, (m,n), α, β, γ, δ, (1,2)) plan_th_jac2jac!(arr::AbstractArray{T}, α, β, γ, δ, dims...) where T = plan_th_jac2jac!(T, size(arr), α, β, γ, δ, dims...) -th_jac2jac(v, α, β, γ, δ, dims...) = plan_th_jac2jac!(eltype(v), size(v), α, β, γ, δ, dims...)*copy(v) \ No newline at end of file +th_jac2jac(v, α, β, γ, δ, dims...) = plan_th_jac2jac!(eltype(v), size(v), α, β, γ, δ, dims...)*copy(v) + + +#### +# cheb2jac +#### + +struct Cheb2JacPlanTH{T, Pl<:Jac2JacPlanTH{T}} <: Plan{T} + jac2jac::Pl +end + + +struct Jac2ChebPlanTH{T, Pl<:Jac2JacPlanTH{T}} <: Plan{T} + jac2jac::Pl +end + + +function jac_cheb_recurrencecoefficients(T, N) + n = 0:N + h = one(T)/2 + A = (2n .+ one(T)) ./ (n .+ one(T)) + A[1] /= 2 + A, Zeros(n), + ((n .- h) .* (n .- h) .* (2n .+ one(T))) ./ ((n .+ one(T)) .* n .* (2n .- one(T))) +end + + +function *(P::Cheb2JacPlanTH{T}, X::AbstractArray) where T + A,B,C = jac_cheb_recurrencecoefficients(T, max(size(X)...)) + + for d in P.jac2jac.dims + if d == 1 + p = forwardrecurrence(size(X,1), A,B,C, one(T)) + X .= p .\ X + else + @assert d == 2 + n = size(X,2) + p = forwardrecurrence(size(X,2), A,B,C, one(T)) + X .= X ./ transpose(p) + end + end + P.jac2jac*X +end + +function *(P::Jac2ChebPlanTH{T}, X::AbstractArray) where T + X = P.jac2jac*X + A,B,C = jac_cheb_recurrencecoefficients(T, max(size(X)...)) + + for d in P.jac2jac.dims + if d == 1 + p = forwardrecurrence(size(X,1), A,B,C, one(T)) + X .= p .* X + else + @assert d == 2 + n = size(X,2) + p = forwardrecurrence(size(X,2), A,B,C, one(T)) + X .= X .* transpose(p) + end + end + X +end + +plan_th_cheb2jac!(::Type{T}, mn, α, β, dims...) where T = Cheb2JacPlanTH(plan_th_jac2jac!(T, mn, -one(α)/2, -one(α)/2, α, β, dims...)) +plan_th_cheb2jac!(arr::AbstractArray{T}, α, β, dims...) where T = plan_th_cheb2jac!(T, size(arr), α, β, dims...) +th_cheb2jac(v, α, β, dims...) = plan_th_cheb2jac!(eltype(v), size(v), α, β, dims...)*copy(v) + +plan_th_jac2cheb!(::Type{T}, mn, α, β, dims...) where T = Jac2ChebPlanTH(plan_th_jac2jac!(T, mn, α, β, -one(α)/2, -one(α)/2, dims...)) +plan_th_jac2cheb!(arr::AbstractArray{T}, α, β, dims...) where T = plan_th_jac2cheb!(T, size(arr), α, β, dims...) +th_jac2cheb(v, α, β, dims...) = plan_th_jac2cheb!(eltype(v), size(v), α, β, dims...)*copy(v) \ No newline at end of file diff --git a/test/toeplitzhankeltests.jl b/test/toeplitzhankeltests.jl index ce1d2f69..b72232ed 100644 --- a/test/toeplitzhankeltests.jl +++ b/test/toeplitzhankeltests.jl @@ -1,7 +1,10 @@ -using FastTransforms, Test +using FastTransforms, Test, Random import FastTransforms: th_leg2cheb, th_cheb2leg, th_leg2chebu, th_ultra2ultra,th_jac2jac, th_leg2chebu, lib_leg2cheb, lib_cheb2leg, lib_ultra2ultra, lib_jac2jac, - plan_th_cheb2leg!, plan_th_leg2chebu!, plan_th_leg2cheb!, plan_th_ultra2ultra!, plan_th_jac2jac! + plan_th_cheb2leg!, plan_th_leg2chebu!, plan_th_leg2cheb!, plan_th_ultra2ultra!, plan_th_jac2jac!, + th_cheb2jac, th_jac2cheb + +Random.seed!(0) @testset "ToeplitzHankel" begin for x in ([1.0], [1.0,2,3,4,5], [1.0+im,2-3im,3+4im,4-5im,5+10im], collect(1.0:1000)) @@ -26,11 +29,14 @@ import FastTransforms: th_leg2cheb, th_cheb2leg, th_leg2chebu, th_ultra2ultra,th @test th_jac2jac(x, -1/2,-1/2,1/2,0) ≈ lib_jac2jac(x, -1/2,-1/2,1/2,0) @test th_jac2jac(x, -1/2,-1/2,0,1/2) ≈ lib_jac2jac(x, -1/2,-1/2,0,1/2) @test th_jac2jac(x, -3/4,-3/4,0,3/4) ≈ lib_jac2jac(x, -3/4,-3/4,0,3/4) - @test th_jac2jac(x,0, 0, 5, 5) ≈ lib_jac2jac(x, 0, 0, 5, 5) if length(x) < 10 + @test th_jac2jac(x,0, 0, 5, 5) ≈ lib_jac2jac(x, 0, 0, 5, 5) @test th_jac2jac(x, 5, 5, 0, 0) ≈ lib_jac2jac(x, 5, 5, 0, 0) end + @test th_cheb2jac(x, 0.2, 0.3) ≈ cheb2jac(x, 0.2, 0.3) + @test th_jac2cheb(x, 0.2, 0.3) ≈ jac2cheb(x, 0.2, 0.3) + @test th_cheb2leg(th_leg2cheb(x)) ≈ x @test th_leg2cheb(th_cheb2leg(x)) ≈ x @test th_ultra2ultra(th_ultra2ultra(x, 0.1, 0.6), 0.6, 0.1) ≈ x @@ -93,6 +99,18 @@ import FastTransforms: th_leg2cheb, th_cheb2leg, th_leg2chebu, th_ultra2ultra,th @test th_jac2jac(X, 0.1, 0.6, 3.1, 2.8, 1) ≈ hcat([jac2jac(X[:,j], 0.1, 0.6, 3.1, 2.8) for j=1:size(X,2)]...) @test th_jac2jac(X, 0.1, 0.6, 3.1, 2.8, 2) ≈ vcat([permutedims(jac2jac(X[k,:], 0.1, 0.6, 3.1, 2.8)) for k=1:size(X,1)]...) @test th_jac2jac(X, 0.1, 0.6, 3.1, 2.8) ≈ th_jac2jac(th_jac2jac(X, 0.1, 0.6, 3.1, 2.8, 1), 0.1, 0.6, 3.1, 2.8, 2) + + @test th_jac2jac(X, -0.5, -0.5, 3.1, 2.8, 1) ≈ hcat([jac2jac(X[:,j], -0.5, -0.5, 3.1, 2.8) for j=1:size(X,2)]...) + @test th_jac2jac(X, -0.5, -0.5, 3.1, 2.8, 2) ≈ vcat([permutedims(jac2jac(X[k,:], -0.5, -0.5, 3.1, 2.8)) for k=1:size(X,1)]...) + @test th_jac2jac(X, -0.5, -0.5, 3.1, 2.8) ≈ th_jac2jac(th_jac2jac(X, -0.5, -0.5, 3.1, 2.8, 1), -0.5, -0.5, 3.1, 2.8, 2) + + @test th_cheb2jac(X, 3.1, 2.8, 1) ≈ hcat([cheb2jac(X[:,j], 3.1, 2.8) for j=1:size(X,2)]...) + @test th_cheb2jac(X, 3.1, 2.8, 2) ≈ vcat([permutedims(cheb2jac(X[k,:], 3.1, 2.8)) for k=1:size(X,1)]...) + @test th_cheb2jac(X, 3.1, 2.8) ≈ th_cheb2jac(th_cheb2jac(X, 3.1, 2.8, 1), 3.1, 2.8, 2) + + @test th_jac2cheb(X, 3.1, 2.8, 1) ≈ hcat([jac2cheb(X[:,j], 3.1, 2.8) for j=1:size(X,2)]...) + @test th_jac2cheb(X, 3.1, 2.8, 2) ≈ vcat([permutedims(jac2cheb(X[k,:], 3.1, 2.8)) for k=1:size(X,1)]...) + @test th_jac2cheb(X, 3.1, 2.8) ≈ th_jac2cheb(th_jac2cheb(X, 3.1, 2.8, 1), 3.1, 2.8, 2) end @testset "BigFloat" begin From 649f25eefbd330390a52b93be4ae61bfcd0156dd Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 22 Oct 2023 09:25:09 +0100 Subject: [PATCH 05/25] Fix bugs hit by th_cheb2jac(x, 1,1) (#230) * Fix bugs hit by th_cheb2jac(x, 1,1) * Update toeplitzhankel.jl --- Project.toml | 2 +- src/toeplitzhankel.jl | 25 +++++++++++++++++++------ test/toeplitzhankeltests.jl | 2 ++ 3 files changed, 22 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index 3a97bbd3..96fdc0a6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.15.9" +version = "0.15.10" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/src/toeplitzhankel.jl b/src/toeplitzhankel.jl index 3bd09d55..4e0ab3db 100644 --- a/src/toeplitzhankel.jl +++ b/src/toeplitzhankel.jl @@ -283,8 +283,17 @@ end isapproxinteger(::Integer) = true isapproxinteger(x) = isinteger(x) || x ≈ round(Int,x) || x+1 ≈ round(Int,x+1) +""" + _nearest_jacobi_par(α, γ) -_nearest_jacobi_par(α, γ) = isapproxinteger(α-γ) ? α : round(Int,α,RoundDown) + mod(γ,1) +returns a number that is an integer different than γ but less than 1 away from α. +""" +function _nearest_jacobi_par(α::T, γ::T) where T + ret = isapproxinteger(α-γ) ? α : round(Int,α,RoundDown) + mod(γ,1) + ret ≤ -1 ? ret + 1 : ret +end +_nearest_jacobi_par(α::T, ::T) where T<:Integer = α +_nearest_jacobi_par(α, γ) = _nearest_jacobi_par(promote(α,γ)...) struct Ultra2UltraPlanTH{T, Plans, Dims} <: Plan{T} @@ -296,12 +305,16 @@ end function *(P::Ultra2UltraPlanTH, A::AbstractArray) ret = A - for p in P.plans - ret = p*ret - end - c = _nearest_jacobi_par(P.λ₁, P.λ₂) + if isapproxinteger(P.λ₂ - P.λ₁) + _ultra2ultra_integerinc!(ret, P.λ₁, P.λ₂, P.dims) + else + for p in P.plans + ret = p*ret + end + c = _nearest_jacobi_par(P.λ₁, P.λ₂) - _ultra2ultra_integerinc!(ret, c, P.λ₂, P.dims) + _ultra2ultra_integerinc!(ret, c, P.λ₂, P.dims) + end end function _ultra2ultraTH_TLC(::Type{S}, mn, λ₁, λ₂, d) where {S} diff --git a/test/toeplitzhankeltests.jl b/test/toeplitzhankeltests.jl index b72232ed..19274747 100644 --- a/test/toeplitzhankeltests.jl +++ b/test/toeplitzhankeltests.jl @@ -36,6 +36,8 @@ Random.seed!(0) @test th_cheb2jac(x, 0.2, 0.3) ≈ cheb2jac(x, 0.2, 0.3) @test th_jac2cheb(x, 0.2, 0.3) ≈ jac2cheb(x, 0.2, 0.3) + @test th_cheb2jac(x, 1, 1) ≈ cheb2jac(x, 1, 1) + @test th_jac2cheb(x, 1, 1) ≈ jac2cheb(x, 1, 1) @test th_cheb2leg(th_leg2cheb(x)) ≈ x @test th_leg2cheb(th_cheb2leg(x)) ≈ x From d9b615df1be09c73dc8d52e73aed73574da61369 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 22 Oct 2023 19:23:14 +0100 Subject: [PATCH 06/25] Matrix Chebyshev U Transforms (#231) * Matrix Chebyshev U Transforms * Cheb U 2 changes * fix inferrability * tests pass * matrix ichebyshevu * tests pass --- Project.toml | 2 +- src/chebyshevtransform.jl | 245 ++++++++++++++++++++++++++------------ test/chebyshevtests.jl | 28 ++++- 3 files changed, 195 insertions(+), 80 deletions(-) diff --git a/Project.toml b/Project.toml index 96fdc0a6..3de40b67 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.15.10" +version = "0.15.11" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/src/chebyshevtransform.jl b/src/chebyshevtransform.jl index 11e07010..28d77eb7 100644 --- a/src/chebyshevtransform.jl +++ b/src/chebyshevtransform.jl @@ -19,9 +19,12 @@ ChebyshevTransformPlan{T,kind}(plan::FFTW.r2rFFTWPlan{T,K,inplace,N,R}) where {T ChebyshevTransformPlan{T,kind,K,inplace,N,R}(plan) # jump through some hoops to make inferrable + +_fftKtype(::Val{N}, _...) where N = Vector{Int32} + function plan_chebyshevtransform!(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - ChebyshevTransformPlan{T,1,Vector{Int32},true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + ChebyshevTransformPlan{T,1,_fftKtype(Val{N}(), dims...),true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else ChebyshevTransformPlan{T,1}(FFTW.plan_r2r!(x, FIRSTKIND, dims...; kws...)) end @@ -34,7 +37,7 @@ end function plan_chebyshevtransform(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - ChebyshevTransformPlan{T,1,Vector{Int32},false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + ChebyshevTransformPlan{T,1,_fftKtype(Val{N}(), dims...),false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else ChebyshevTransformPlan{T,1}(FFTW.plan_r2r(x, FIRSTKIND, dims...; kws...)) end @@ -149,17 +152,15 @@ end end function *(P::ChebyshevTransformPlan{T,1,K,true,N}, x::AbstractArray{T,N}) where {T,K,N} - n = length(x) - n == 0 && return x + isempty(x) && return x y = P.plan*x # will be === x if in-place _cheb1_rescale!(P.plan.region, y) end function mul!(y::AbstractArray{T,N}, P::ChebyshevTransformPlan{T,1,K,false,N}, x::AbstractArray{<:Any,N}) where {T,K,N} - n = length(x) - length(y) == n || throw(DimensionMismatch("output must match dimension")) - n == 0 && return y + size(y) == size(x) || throw(DimensionMismatch("output must match dimension")) + isempty(x) && return y _plan_mul!(y, P.plan, x) _cheb1_rescale!(P.plan.region, y) end @@ -248,7 +249,7 @@ inv(P::IChebyshevTransformPlan{T,1}) where {T} = ChebyshevTransformPlan{T,1}(inv function plan_ichebyshevtransform!(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - IChebyshevTransformPlan{T,1,Vector{Int32},true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + IChebyshevTransformPlan{T,1,_fftKtype(Val{N}(), dims...),true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else IChebyshevTransformPlan{T,1}(FFTW.plan_r2r!(x, IFIRSTKIND, dims...; kws...)) end @@ -260,7 +261,7 @@ end function plan_ichebyshevtransform(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - IChebyshevTransformPlan{T,1,Vector{Int32},false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + IChebyshevTransformPlan{T,1,_fftKtype(Val{N}(), dims...),false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else IChebyshevTransformPlan{T,1}(FFTW.plan_r2r(x, IFIRSTKIND, dims...; kws...)) end @@ -308,7 +309,7 @@ function mul!(y::AbstractArray{T,N}, P::IChebyshevTransformPlan{T,1,K,false,N}, size(y) == size(x) || throw(DimensionMismatch("output must match dimension")) isempty(x) && return y - _icheb1_prescale!(P.plan.region, x) # Todo: don't mutate x + _icheb1_prescale!(P.plan.region, x) # TODO: don't mutate x _plan_mul!(y, P.plan, x) _icheb1_postscale!(P.plan.region, x) ldiv!(2^length(P.plan.region), y) @@ -371,7 +372,9 @@ ichebyshevtransform!(x::AbstractArray, dims...; kwds...) = plan_ichebyshevtransf ichebyshevtransform(x, dims...; kwds...) = plan_ichebyshevtransform(x, dims...; kwds...)*x -## Chebyshev U +####### +# Chebyshev U +####### const UFIRSTKIND = FFTW.RODFT10 const USECONDKIND = FFTW.RODFT00 @@ -388,7 +391,7 @@ ChebyshevUTransformPlan{T,kind}(plan::FFTW.r2rFFTWPlan{T,K,inplace,N,R}) where { function plan_chebyshevutransform!(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - ChebyshevUTransformPlan{T,1,Vector{Int32},true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + ChebyshevUTransformPlan{T,1,_fftKtype(Val{N}(), dims...),true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else ChebyshevUTransformPlan{T,1}(FFTW.plan_r2r!(x, UFIRSTKIND, dims...; kws...)) end @@ -400,13 +403,19 @@ end function plan_chebyshevutransform(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - ChebyshevUTransformPlan{T,1,Vector{Int32},false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + ChebyshevUTransformPlan{T,1,_fftKtype(Val{N}(), dims...),false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else ChebyshevUTransformPlan{T,1}(FFTW.plan_r2r(x, UFIRSTKIND, dims...; kws...)) end end function plan_chebyshevutransform(x::AbstractArray{T,N}, ::Val{2}, dims...; kws...) where {T<:fftwNumber,N} - any(≤(1),size(x)) && throw(ArgumentError("Array must contain at least 2 entries")) + if isempty(dims) + any(≤(1), size(x)) && throw(ArgumentError("Array must contain at least 2 entries")) + else + for d in dims[1] + size(x,d) ≤ 1 && throw(ArgumentError("Array must contain at least 2 entries")) + end + end ChebyshevUTransformPlan{T,2}(FFTW.plan_r2r(x, USECONDKIND, dims...; kws...)) end @@ -414,75 +423,141 @@ plan_chebyshevutransform!(x::AbstractArray, dims...; kws...) = plan_chebyshevutr plan_chebyshevutransform(x::AbstractArray, dims...; kws...) = plan_chebyshevutransform(x, Val(1), dims...; kws...) -@inline function _chebu1_prescale!(_, x::AbstractVector{T}) where T - n = length(x) - for k=1:n # sqrt(1-x_j^2) weight - x[k] *= sinpi(one(T)/(2n) + (k-one(T))/n)/n +@inline function _chebu1_prescale!(d::Number, x::AbstractVecOrMat{T}) where T + m,n = size(x,1),size(x,2) + if d == 1 + for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight + x[k,j] *= sinpi(one(T)/(2m) + (k-one(T))/m)/m + end + else + @assert d == 2 + for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight + x[k,j] *= sinpi(one(T)/(2n) + (j-one(T))/n)/n + end end x end -@inline function _chebu1_postscale!(_, x::AbstractVector{T}) where T - n = length(x) - for k=1:n # sqrt(1-x_j^2) weight - x[k] /= sinpi(one(T)/(2n) + (k-one(T))/n)/n +@inline function _chebu1_prescale!(d, y::AbstractArray) + for k in d + _chebu1_prescale!(k, y) + end + y +end + +@inline function _chebu1_postscale!(d::Number, x::AbstractVecOrMat{T}) where T + m,n = size(x,1),size(x,2) + if d == 1 + for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight + x[k,j] /= sinpi(one(T)/(2m) + (k-one(T))/m)/m + end + else + @assert d == 2 + for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight + x[k,j] /= sinpi(one(T)/(2n) + (j-one(T))/n)/n + end end x end -function *(P::ChebyshevUTransformPlan{T,1,K,true}, x::AbstractVector{T}) where {T,K} +@inline function _chebu1_postscale!(d, y::AbstractArray) + for k in d + _chebu1_postscale!(k, y) + end + y +end + +function *(P::ChebyshevUTransformPlan{T,1,K,true,N}, x::AbstractArray{T,N}) where {T,K,N} length(x) ≤ 1 && return x _chebu1_prescale!(P.plan.region, x) P.plan * x end -function mul!(y::AbstractVector{T}, P::ChebyshevUTransformPlan{T,1,K,false}, x::AbstractVector{T}) where {T,K} - n = length(x) - length(x) ≤ 1 && return copyto!(y, x) - _chebu1_prescale!(P.plan.region, x) +function mul!(y::AbstractArray{T}, P::ChebyshevUTransformPlan{T,1,K,false}, x::AbstractArray{T}) where {T,K} + size(y) == size(x) || throw(DimensionMismatch("output must match dimension")) + isempty(x) && return y + _chebu1_prescale!(P.plan.region, x) # Todo don't mutate x _plan_mul!(y, P.plan, x) _chebu1_postscale!(P.plan.region, x) + for d in P.plan.region + size(y,d) == 1 && ldiv!(2, y) # fix doubling + end y end -@inline function _chebu2_prescale!(_, x::AbstractVector{T}) where T - n = length(x) - c = one(T)/ (n+1) - for k=1:n # sqrt(1-x_j^2) weight - x[k] *= sinpi(k*c) +@inline function _chebu2_prescale!(d::Number, x::AbstractVecOrMat{T}) where T + m,n = size(x,1),size(x,2) + if d == 1 + c = one(T)/ (m+1) + for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight + x[k,j] *= sinpi(k*c) + end + else + @assert d == 2 + c = one(T)/ (n+1) + for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight + x[k,j] *= sinpi(j*c) + end end x end -@inline function _chebu2_postscale!(_, x::AbstractVector{T}) where T - n = length(x) - c = one(T)/ (n+1) - @inbounds for k=1:n # sqrt(1-x_j^2) weight - x[k] /= sinpi(k*c) +@inline function _chebu2_prescale!(d, y::AbstractArray) + for k in d + _chebu2_prescale!(k, y) + end + y +end + + +@inline function _chebu2_postscale!(d::Number, x::AbstractVecOrMat{T}) where T + m,n = size(x,1),size(x,2) + if d == 1 + c = one(T)/ (m+1) + for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight + x[k,j] /= sinpi(k*c) + end + else + @assert d == 2 + c = one(T)/ (n+1) + for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight + x[k,j] /= sinpi(j*c) + end end x end -function *(P::ChebyshevUTransformPlan{T,2,K,true}, x::AbstractVector{T}) where {T,K} - n = length(x) - n ≤ 1 && return x - _chebu2_prescale!(P.plan.region, x) - lmul!(one(T)/ (n+1), P.plan * x) +@inline function _chebu2_postscale!(d, y::AbstractArray) + for k in d + _chebu2_postscale!(k, y) + end + y end -function mul!(y::AbstractVector{T}, P::ChebyshevUTransformPlan{T,2,K,false}, x::AbstractVector{T}) where {T,K} - n = length(x) - n ≤ 1 && return copyto!(y, x) +function *(P::ChebyshevUTransformPlan{T,2,K,true,N}, x::AbstractArray{T,N}) where {T,K,N} + sc = one(T) + for d in P.plan.region + sc *= one(T)/(size(x,d)+1) + end _chebu2_prescale!(P.plan.region, x) + lmul!(sc, P.plan * x) +end + +function mul!(y::AbstractArray{T}, P::ChebyshevUTransformPlan{T,2,K,false}, x::AbstractArray{T}) where {T,K} + sc = one(T) + for d in P.plan.region + sc *= one(T)/(size(x,d)+1) + end + _chebu2_prescale!(P.plan.region, x) # TODO don't mutate x _plan_mul!(y, P.plan, x) _chebu2_postscale!(P.plan.region, x) - lmul!(one(T)/ (n+1), y) + lmul!(sc, y) end *(P::ChebyshevUTransformPlan{T,kind,K,false,N}, x::AbstractArray{T,N}) where {T,kind,K,N} = mul!(similar(x), P, x) -chebyshevutransform!(x::AbstractVector{T}, dims...; kws...) where {T<:fftwNumber} = +chebyshevutransform!(x::AbstractArray{T}, dims...; kws...) where {T<:fftwNumber} = plan_chebyshevutransform!(x, dims...; kws...)*x @@ -509,19 +584,19 @@ IChebyshevUTransformPlan{T,kind}(F::FFTW.r2rFFTWPlan{T,K,inplace,N,R}) where {T, function plan_ichebyshevutransform!(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - IChebyshevUTransformPlan{T,1,Vector{Int32},true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + IChebyshevUTransformPlan{T,1,_fftKtype(Val{N}(), dims...),true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else IChebyshevUTransformPlan{T,1}(FFTW.plan_r2r!(x, IUFIRSTKIND, dims...; kws...)) end end function plan_ichebyshevutransform!(x::AbstractArray{T,N}, ::Val{2}, dims...; kws...) where {T<:fftwNumber,N} any(≤(1),size(x)) && throw(ArgumentError("Array must contain at least 2 entries")) - IChebyshevUTransformPlan{T,2}(FFTW.plan_r2r!(x, USECONDKIND)) + IChebyshevUTransformPlan{T,2}(FFTW.plan_r2r!(x, USECONDKIND, dims...)) end function plan_ichebyshevutransform(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - IChebyshevUTransformPlan{T,1,Vector{Int32},false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + IChebyshevUTransformPlan{T,1,_fftKtype(Val{N}(), dims...),false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else IChebyshevUTransformPlan{T,1}(FFTW.plan_r2r(x, IUFIRSTKIND, dims...; kws...)) end @@ -543,42 +618,59 @@ inv(P::IChebyshevUTransformPlan{T,2}) where {T} = ChebyshevUTransformPlan{T,2}(P inv(P::ChebyshevUTransformPlan{T,1}) where {T} = IChebyshevUTransformPlan{T,1}(inv(P.plan).p) inv(P::IChebyshevUTransformPlan{T,1}) where {T} = ChebyshevUTransformPlan{T,1}(inv(P.plan).p) - -function _ichebyu1_postscale!(_, x::AbstractVector{T}) where T - n = length(x) - @inbounds for k=1:n # sqrt(1-x_j^2) weight - x[k] /= 2sinpi(one(T)/(2n) + (k-one(T))/n) +@inline function _ichebu1_postscale!(d::Number, x::AbstractVecOrMat{T}) where T + m,n = size(x,1),size(x,2) + if d == 1 + for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight + x[k,j] /= 2sinpi(one(T)/(2m) + (k-one(T))/m) + end + else + @assert d == 2 + for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight + x[k,j] /= 2sinpi(one(T)/(2n) + (j-one(T))/n) + end end x end -function *(P::IChebyshevUTransformPlan{T,1,K,true}, x::AbstractVector{T}) where {T<:fftwNumber,K} - n = length(x) - n ≤ 1 && return x - x = P.plan * x - _ichebyu1_postscale!(P.plan.region, x) +@inline function _ichebu1_postscale!(d, y::AbstractArray) + for k in d + _ichebu1_postscale!(k, y) + end + y end -function mul!(y::AbstractVector{T}, P::IChebyshevUTransformPlan{T,1,K,false}, x::AbstractVector{T}) where {T<:fftwNumber,K} - n = length(x) - length(y) == n || throw(DimensionMismatch("output must match dimension")) - n ≤ 1 && return x +function *(P::IChebyshevUTransformPlan{T,1,K,true}, x::AbstractArray{T}) where {T<:fftwNumber,K} + length(x) ≤ 1 && return x + x = P.plan * x + _ichebu1_postscale!(P.plan.region, x) +end +function mul!(y::AbstractArray{T}, P::IChebyshevUTransformPlan{T,1,K,false}, x::AbstractArray{T}) where {T<:fftwNumber,K} + size(y) == size(x) || throw(DimensionMismatch("output must match dimension")) + isempty(x) && return y _plan_mul!(y, P.plan, x) - _ichebyu1_postscale!(P.plan.region, y) + _ichebu1_postscale!(P.plan.region, y) + for d in P.plan.region + size(y,d) == 1 && lmul!(2, y) # fix doubling + end + y end -function _ichebu2_rescale!(_, x::AbstractVector{T}) where T - n = length(x) - c = one(T)/ (n+1) - for k=1:n # sqrt(1-x_j^2) weight - x[k] /= sinpi(k*c) - end +function _ichebu2_rescale!(d::Number, x::AbstractArray{T}) where T + _chebu2_postscale!(d, x) ldiv!(2, x) x end -function *(P::IChebyshevUTransformPlan{T,2,K,true}, x::AbstractVector{T}) where {T<:fftwNumber,K} +@inline function _ichebu2_rescale!(d, y::AbstractArray) + for k in d + _ichebu2_rescale!(k, y) + end + y +end + +function *(P::IChebyshevUTransformPlan{T,2,K,true}, x::AbstractArray{T}) where {T<:fftwNumber,K} n = length(x) n ≤ 1 && return x @@ -586,16 +678,15 @@ function *(P::IChebyshevUTransformPlan{T,2,K,true}, x::AbstractVector{T}) where _ichebu2_rescale!(P.plan.region, x) end -function mul!(y::AbstractVector{T}, P::IChebyshevUTransformPlan{T,2,K,false}, x::AbstractVector{T}) where {T<:fftwNumber,K} - n = length(x) - length(y) == n || throw(DimensionMismatch("output must match dimension")) - n ≤ 1 && return x +function mul!(y::AbstractArray{T}, P::IChebyshevUTransformPlan{T,2,K,false}, x::AbstractArray{T}) where {T<:fftwNumber,K} + size(y) == size(x) || throw(DimensionMismatch("output must match dimension")) + length(x) ≤ 1 && return x _plan_mul!(y, P.plan, x) _ichebu2_rescale!(P.plan.region, y) end -ichebyshevutransform!(x::AbstractVector{T}, dims...; kwds...) where {T<:fftwNumber} = +ichebyshevutransform!(x::AbstractArray{T}, dims...; kwds...) where {T<:fftwNumber} = plan_ichebyshevutransform!(x, dims...; kwds...)*x ichebyshevutransform(x, dims...; kwds...) = plan_ichebyshevutransform(x, dims...; kwds...)*x diff --git a/test/chebyshevtests.jl b/test/chebyshevtests.jl index 696c13f0..82b8da5c 100644 --- a/test/chebyshevtests.jl +++ b/test/chebyshevtests.jl @@ -154,6 +154,7 @@ using FastTransforms, Test p_1 = chebyshevpoints(T, n) f = exp.(p_1) g = @inferred(chebyshevutransform(f)) + @test f ≈ exp.(p_1) f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-1]' * g @test f̃(0.1) ≈ exp(T(0.1)) @@ -221,7 +222,7 @@ using FastTransforms, Test f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-3]' * g @test f̃(0.1) ≈ exp(T(0.1)) - @test @inferred(ichebyshevutransform(g, Val(2))) ≈ exp.(p_2) + @test @inferred(ichebyshevutransform(g, Val(2))) ≈ f ≈ exp.(p_2) fcopy = copy(f) gcopy = copy(g) @@ -291,6 +292,29 @@ using FastTransforms, Test @test chebyshevtransform(ichebyshevtransform(X)) ≈ X end + @testset "chebyshevutransform" begin + @test @inferred(chebyshevutransform(X,1)) ≈ @inferred(chebyshevutransform!(copy(X),1)) ≈ hcat(chebyshevutransform.([X[:,k] for k=axes(X,2)])...) + @test chebyshevutransform(X,2) ≈ chebyshevutransform!(copy(X),2) ≈ hcat(chebyshevutransform.([X[k,:] for k=axes(X,1)])...)' + @test @inferred(chebyshevutransform(X,Val(2),1)) ≈ @inferred(chebyshevutransform!(copy(X),Val(2),1)) ≈ hcat(chebyshevutransform.([X[:,k] for k=axes(X,2)],Val(2))...) + @test chebyshevutransform(X,Val(2),2) ≈ chebyshevutransform!(copy(X),Val(2),2) ≈ hcat(chebyshevutransform.([X[k,:] for k=axes(X,1)],Val(2))...)' + + @test @inferred(chebyshevutransform(X)) ≈ @inferred(chebyshevutransform!(copy(X))) ≈ chebyshevutransform(chebyshevutransform(X,1),2) + @test @inferred(chebyshevutransform(X,Val(2))) ≈ @inferred(chebyshevutransform!(copy(X),Val(2))) ≈ chebyshevutransform(chebyshevutransform(X,Val(2),1),Val(2),2) + end + + @testset "ichebyshevutransform" begin + @test @inferred(ichebyshevutransform(X,1)) ≈ @inferred(ichebyshevutransform!(copy(X),1)) ≈ hcat(ichebyshevutransform.([X[:,k] for k=axes(X,2)])...) + @test ichebyshevutransform(X,2) ≈ ichebyshevutransform!(copy(X),2) ≈ hcat(ichebyshevutransform.([X[k,:] for k=axes(X,1)])...)' + @test @inferred(ichebyshevutransform(X,Val(2),1)) ≈ @inferred(ichebyshevutransform!(copy(X),Val(2),1)) ≈ hcat(ichebyshevutransform.([X[:,k] for k=axes(X,2)],Val(2))...) + @test ichebyshevutransform(X,Val(2),2) ≈ ichebyshevutransform!(copy(X),Val(2),2) ≈ hcat(ichebyshevutransform.([X[k,:] for k=axes(X,1)],Val(2))...)' + + @test @inferred(ichebyshevutransform(X)) ≈ @inferred(ichebyshevutransform!(copy(X))) ≈ ichebyshevutransform(ichebyshevutransform(X,1),2) + @test @inferred(ichebyshevutransform(X,Val(2))) ≈ @inferred(ichebyshevutransform!(copy(X),Val(2))) ≈ ichebyshevutransform(ichebyshevutransform(X,Val(2),1),Val(2),2) + + @test ichebyshevutransform(chebyshevutransform(X)) ≈ X + @test chebyshevutransform(ichebyshevutransform(X)) ≈ X + end + X = randn(1,1) @test chebyshevtransform!(copy(X), Val(1)) == ichebyshevtransform!(copy(X), Val(1)) == X @test_throws ArgumentError chebyshevtransform!(copy(X), Val(2)) @@ -394,7 +418,7 @@ using FastTransforms, Test plan_chebyshevutransform(X,Val(1),2), plan_chebyshevutransform(X, Val(2),2), plan_ichebyshevutransform(X,Val(1),1), plan_ichebyshevutransform(X, Val(2),1), plan_ichebyshevutransform(X,Val(1),2), plan_ichebyshevutransform(X, Val(2),2)) - @test_broken F \ (F*X) ≈ F * (F\X) ≈ X + @test F \ (F*X) ≈ F * (F\X) ≈ X end end end From db873136f7b2eaf5f8ee33701359ad49dc7443d5 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 7 Nov 2023 16:21:31 +0530 Subject: [PATCH 07/25] Docstrings for orthogonal polynomial transforms (#199) * Docstrings for orthogonal polynomial transforms * Update leg2cheb Co-authored-by: Sheehan Olver * Update cheb2leg * StridedVector -> AbstractVector --------- Co-authored-by: Sheehan Olver --- src/FastTransforms.jl | 1 + src/docstrings.jl | 121 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 122 insertions(+) create mode 100644 src/docstrings.jl diff --git a/src/FastTransforms.jl b/src/FastTransforms.jl index c8017b27..f9d99e74 100644 --- a/src/FastTransforms.jl +++ b/src/FastTransforms.jl @@ -134,5 +134,6 @@ end # end # end +include("docstrings.jl") end # module diff --git a/src/docstrings.jl b/src/docstrings.jl new file mode 100644 index 00000000..c3ecd7a3 --- /dev/null +++ b/src/docstrings.jl @@ -0,0 +1,121 @@ +""" + leg2cheb(v::AbstractVector; normleg::Bool=false, normcheb::Bool=false) + +Convert the vector of expansions coefficients `v` from a Legendre to a Chebyshev basis. +The keyword arguments denote whether the bases are normalized. +""" +leg2cheb + +""" + cheb2leg(v::AbstractVector; normcheb::Bool=false, normleg::Bool=false) + +Convert the vector of expansions coefficients `v` from a Chebyshev to a Legendre basis. +The keyword arguments denote whether the bases are normalized. +""" +cheb2leg + +""" + ultra2ultra(v::AbstractVector, λ, μ; norm1::Bool=false, norm2::Bool=false) + +Convert the vector of expansions coefficients `v` from an Ultraspherical basis of +order `λ` to an Ultraspherical basis of order `μ`. +The keyword arguments denote whether the bases are normalized. +""" +ultra2ultra + +""" + jac2jac(v::AbstractVector, α, β, γ, δ; norm1::Bool=false, norm2::Bool=false) + +Convert the vector of expansions coefficients `v` from a Jacobi basis of +order `(α,β)` to a Jacobi basis of order `(γ,δ)`. +The keyword arguments denote whether the bases are normalized. +""" +jac2jac + +""" + lag2lag(v::AbstractVector, α, β; norm1::Bool=false, norm2::Bool=false) + +Convert the vector of expansions coefficients `v` from a Laguerre basis of +order `α` to a La basis of order `β`. +The keyword arguments denote whether the bases are normalized.""" +lag2lag + +""" + jac2ultra(v::AbstractVector, α, β, λ; normjac::Bool=false, normultra::Bool=false) + +Convert the vector of expansions coefficients `v` from a Jacobi basis of +order `(α,β)` to an Ultraspherical basis of order `λ`. +The keyword arguments denote whether the bases are normalized.""" +jac2ultra + +""" + ultra2jac(v::AbstractVector, λ, α, β; normultra::Bool=false, normjac::Bool=false) + +Convert the vector of expansions coefficients `v` from an Ultraspherical basis of +order `λ` to a Jacobi basis of order `(α,β)`. +The keyword arguments denote whether the bases are normalized. +""" +ultra2jac + +""" + jac2cheb(v::AbstractVector, α, β; normjac::Bool=false, normcheb::Bool=false) + +Convert the vector of expansions coefficients `v` from a Jacobi basis of +order `(α,β)` to a Chebyshev basis. +The keyword arguments denote whether the bases are normalized. +""" +jac2cheb + +""" + cheb2jac(v::AbstractVector, α, β; normcheb::Bool=false, normjac::Bool=false) + +Convert the vector of expansions coefficients `v` from a Chebyshev basis to a +Jacobi basis of order `(α,β)`. +The keyword arguments denote whether the bases are normalized. +""" +cheb2jac + +""" + ultra2cheb(v::AbstractVector, λ; normultra::Bool=false, normcheb::Bool=false) + +Convert the vector of expansions coefficients `v` from an Ultraspherical basis of +order `λ` to a Chebyshev basis. +The keyword arguments denote whether the bases are normalized. +""" +ultra2cheb + +""" + cheb2ultra(v::AbstractVector, λ; normcheb::Bool=false, normultra::Bool=false) + +Convert the vector of expansions coefficients `v` from a Chebyshev basis +to an Ultraspherical basis of order `λ`. +The keyword arguments denote whether the bases are normalized. +""" +cheb2ultra + +""" + associatedjac2jac(v::AbstractVector, c::Integer, α, β, γ, δ; norm1::Bool=false, norm2::Bool=false) + +Convert the vector of expansions coefficients `v` from an associated Jacobi basis +of orders `(α,β)` to a Jacobi basis of order `(γ,δ)`. +The keyword arguments denote whether the bases are normalized. +""" +associatedjac2jac + +""" + modifiedjac2jac(v::AbstractVector{T}, α, β, u::Vector{T}; verbose::Bool=false) where {T} + modifiedjac2jac(v::AbstractVector{T}, α, β, u::Vector{T}, v::Vector{T}; verbose::Bool=false) where {T} +""" +modifiedjac2jac + +""" + modifiedlag2lag(v::AbstractVector{T}, α, u::Vector{T}; verbose::Bool=false) + modifiedlag2lag(v::AbstractVector{T}, α, u::Vector{T}, v::Vector{T}; verbose::Bool=false) where {T} +""" +modifiedlag2lag + +""" + modifiedherm2herm(v::AbstractVector{T}, u::Vector{T}; verbose::Bool=false) + modifiedherm2herm(v::AbstractVector{T}, u::Vector{T}, v::Vector{T}; verbose::Bool=false) where {T} +""" +modifiedherm2herm From 198d22016aeca1d5cb260bedc8d8864612ee829e Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Fri, 10 Nov 2023 22:46:21 +0000 Subject: [PATCH 08/25] Support arbitrary tensors in Chebyshev U (#232) * Support arbitrary tensors * tests pass --- src/chebyshevtransform.jl | 100 ++++++++++++-------------------------- test/chebyshevtests.jl | 45 ++++++++++++++++- 2 files changed, 74 insertions(+), 71 deletions(-) diff --git a/src/chebyshevtransform.jl b/src/chebyshevtransform.jl index 28d77eb7..01dd61af 100644 --- a/src/chebyshevtransform.jl +++ b/src/chebyshevtransform.jl @@ -423,19 +423,13 @@ plan_chebyshevutransform!(x::AbstractArray, dims...; kws...) = plan_chebyshevutr plan_chebyshevutransform(x::AbstractArray, dims...; kws...) = plan_chebyshevutransform(x, Val(1), dims...; kws...) -@inline function _chebu1_prescale!(d::Number, x::AbstractVecOrMat{T}) where T - m,n = size(x,1),size(x,2) - if d == 1 - for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight - x[k,j] *= sinpi(one(T)/(2m) + (k-one(T))/m)/m - end - else - @assert d == 2 - for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight - x[k,j] *= sinpi(one(T)/(2n) + (j-one(T))/n)/n - end - end - x +_permfirst(d, N) = [d; 1:d-1; d+1:N] + +@inline function _chebu1_prescale!(d::Number, X::AbstractArray{T,N}) where {T,N} + X̃ = PermutedDimsArray(X, _permfirst(d, N)) + m = size(X̃,1) + X̃ .= (sinpi.(one(T)/(2m) .+ ((1:m) .- one(T))/m) ./ m) .* X̃ + X end @inline function _chebu1_prescale!(d, y::AbstractArray) @@ -445,19 +439,11 @@ end y end -@inline function _chebu1_postscale!(d::Number, x::AbstractVecOrMat{T}) where T - m,n = size(x,1),size(x,2) - if d == 1 - for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight - x[k,j] /= sinpi(one(T)/(2m) + (k-one(T))/m)/m - end - else - @assert d == 2 - for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight - x[k,j] /= sinpi(one(T)/(2n) + (j-one(T))/n)/n - end - end - x +@inline function _chebu1_postscale!(d::Number, X::AbstractArray{T,N}) where {T,N} + X̃ = PermutedDimsArray(X, _permfirst(d, N)) + m = size(X̃,1) + X̃ .= X̃ ./ (sinpi.(one(T)/(2m) .+ ((1:m) .- one(T))/m) ./ m) + X end @inline function _chebu1_postscale!(d, y::AbstractArray) @@ -485,21 +471,13 @@ function mul!(y::AbstractArray{T}, P::ChebyshevUTransformPlan{T,1,K,false}, x::A y end -@inline function _chebu2_prescale!(d::Number, x::AbstractVecOrMat{T}) where T - m,n = size(x,1),size(x,2) - if d == 1 - c = one(T)/ (m+1) - for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight - x[k,j] *= sinpi(k*c) - end - else - @assert d == 2 - c = one(T)/ (n+1) - for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight - x[k,j] *= sinpi(j*c) - end - end - x + +@inline function _chebu2_prescale!(d::Number, X::AbstractArray{T,N}) where {T,N} + X̃ = PermutedDimsArray(X, _permfirst(d, N)) + m = size(X̃,1) + c = one(T)/ (m+1) + X̃ .= sinpi.((1:m) .* c) .* X̃ + X end @inline function _chebu2_prescale!(d, y::AbstractArray) @@ -510,21 +488,12 @@ end end -@inline function _chebu2_postscale!(d::Number, x::AbstractVecOrMat{T}) where T - m,n = size(x,1),size(x,2) - if d == 1 - c = one(T)/ (m+1) - for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight - x[k,j] /= sinpi(k*c) - end - else - @assert d == 2 - c = one(T)/ (n+1) - for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight - x[k,j] /= sinpi(j*c) - end - end - x +@inline function _chebu2_postscale!(d::Number, X::AbstractArray{T,N}) where {T,N} + X̃ = PermutedDimsArray(X, _permfirst(d, N)) + m = size(X̃,1) + c = one(T)/ (m+1) + X̃ .= X̃ ./ sinpi.((1:m) .* c) + X end @inline function _chebu2_postscale!(d, y::AbstractArray) @@ -618,21 +587,14 @@ inv(P::IChebyshevUTransformPlan{T,2}) where {T} = ChebyshevUTransformPlan{T,2}(P inv(P::ChebyshevUTransformPlan{T,1}) where {T} = IChebyshevUTransformPlan{T,1}(inv(P.plan).p) inv(P::IChebyshevUTransformPlan{T,1}) where {T} = ChebyshevUTransformPlan{T,1}(inv(P.plan).p) -@inline function _ichebu1_postscale!(d::Number, x::AbstractVecOrMat{T}) where T - m,n = size(x,1),size(x,2) - if d == 1 - for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight - x[k,j] /= 2sinpi(one(T)/(2m) + (k-one(T))/m) - end - else - @assert d == 2 - for j = 1:n, k = 1:m # sqrt(1-x_j^2) weight - x[k,j] /= 2sinpi(one(T)/(2n) + (j-one(T))/n) - end - end - x +@inline function _ichebu1_postscale!(d::Number, X::AbstractArray{T,N}) where {T,N} + X̃ = PermutedDimsArray(X, _permfirst(d, N)) + m = size(X̃,1) + X̃ .= X̃ ./ (2 .* sinpi.(one(T)/(2m) .+ ((1:m) .- one(T))/m)) + X end + @inline function _ichebu1_postscale!(d, y::AbstractArray) for k in d _ichebu1_postscale!(k, y) diff --git a/test/chebyshevtests.jl b/test/chebyshevtests.jl index 82b8da5c..d6f704bf 100644 --- a/test/chebyshevtests.jl +++ b/test/chebyshevtests.jl @@ -164,11 +164,11 @@ using FastTransforms, Test gcopy = copy(g) P = @inferred(plan_chebyshevutransform(f)) @test P*f ≈ g - @test f == fcopy + @test f ≈ fcopy @test_throws ArgumentError P * T[1,2] P = @inferred(plan_chebyshevutransform(f, 1:1)) @test P*f ≈ g - @test f == fcopy + @test f ≈ fcopy @test_throws ArgumentError P * T[1,2] P = @inferred(plan_chebyshevutransform!(f)) @@ -364,6 +364,47 @@ using FastTransforms, Test @test ichebyshevtransform(chebyshevtransform(X)) ≈ X @test chebyshevtransform(ichebyshevtransform(X)) ≈ X end + + @testset "chebyshevutransform" begin + for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = chebyshevutransform(X[:,k,j]) end + @test @inferred(chebyshevutransform(X,1)) ≈ @inferred(chebyshevutransform!(copy(X),1)) ≈ X̃ + for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = chebyshevutransform(X[k,:,j]) end + @test chebyshevutransform(X,2) ≈ chebyshevutransform!(copy(X),2) ≈ X̃ + for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = chebyshevutransform(X[k,j,:]) end + @test chebyshevutransform(X,3) ≈ chebyshevutransform!(copy(X),3) ≈ X̃ + + for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = chebyshevutransform(X[:,k,j],Val(2)) end + @test @inferred(chebyshevutransform(X,Val(2),1)) ≈ @inferred(chebyshevutransform!(copy(X),Val(2),1)) ≈ X̃ + for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = chebyshevutransform(X[k,:,j],Val(2)) end + @test chebyshevutransform(X,Val(2),2) ≈ chebyshevutransform!(copy(X),Val(2),2) ≈ X̃ + for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = chebyshevutransform(X[k,j,:],Val(2)) end + @test chebyshevutransform(X,Val(2),3) ≈ chebyshevutransform!(copy(X),Val(2),3) ≈ X̃ + + @test @inferred(chebyshevutransform(X)) ≈ @inferred(chebyshevutransform!(copy(X))) ≈ chebyshevutransform(chebyshevutransform(chebyshevutransform(X,1),2),3) + @test @inferred(chebyshevutransform(X,Val(2))) ≈ @inferred(chebyshevutransform!(copy(X),Val(2))) ≈ chebyshevutransform(chebyshevutransform(chebyshevutransform(X,Val(2),1),Val(2),2),Val(2),3) + end + + @testset "ichebyshevutransform" begin + for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = ichebyshevutransform(X[:,k,j]) end + @test @inferred(ichebyshevutransform(X,1)) ≈ @inferred(ichebyshevutransform!(copy(X),1)) ≈ X̃ + for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = ichebyshevutransform(X[k,:,j]) end + @test ichebyshevutransform(X,2) ≈ ichebyshevutransform!(copy(X),2) ≈ X̃ + for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = ichebyshevutransform(X[k,j,:]) end + @test ichebyshevutransform(X,3) ≈ ichebyshevutransform!(copy(X),3) ≈ X̃ + + for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = ichebyshevutransform(X[:,k,j],Val(2)) end + @test @inferred(ichebyshevutransform(X,Val(2),1)) ≈ @inferred(ichebyshevutransform!(copy(X),Val(2),1)) ≈ X̃ + for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = ichebyshevutransform(X[k,:,j],Val(2)) end + @test ichebyshevutransform(X,Val(2),2) ≈ ichebyshevutransform!(copy(X),Val(2),2) ≈ X̃ + for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = ichebyshevutransform(X[k,j,:],Val(2)) end + @test ichebyshevutransform(X,Val(2),3) ≈ ichebyshevutransform!(copy(X),Val(2),3) ≈ X̃ + + @test @inferred(ichebyshevutransform(X)) ≈ @inferred(ichebyshevutransform!(copy(X))) ≈ ichebyshevutransform(ichebyshevutransform(ichebyshevutransform(X,1),2),3) + @test @inferred(ichebyshevutransform(X,Val(2))) ≈ @inferred(ichebyshevutransform!(copy(X),Val(2))) ≈ ichebyshevutransform(ichebyshevutransform(ichebyshevutransform(X,Val(2),1),Val(2),2),Val(2),3) + + @test ichebyshevutransform(chebyshevutransform(X)) ≈ X + @test chebyshevutransform(ichebyshevutransform(X)) ≈ X + end X = randn(1,1,1) @test chebyshevtransform!(copy(X), Val(1)) == ichebyshevtransform!(copy(X), Val(1)) == X From dcaf8f6338a35f0f0eb512dd5d5f323178dbc179 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sun, 12 Nov 2023 21:45:18 +0000 Subject: [PATCH 09/25] Support higher order tensors in ToeplitzHankel (#233) * Support higher order tensors in ToeplitzHankel * Update toeplitzhankel.jl * Generalise ToeplitzPlan * generalis plan_upper for matrices * consolidate plan_uppertoeplitz * 4-tensor tests * Update toeplitzplans.jl * generalise _th_applymul! * Tesnor leg2cheb * tests pass * cheb2leg tensor support * fixes * v0.15.12 * Update toeplitzplans.jl --- Project.toml | 2 +- src/FastTransforms.jl | 2 +- src/toeplitzhankel.jl | 136 +++++++-------------- src/toeplitzplans.jl | 231 ++++++++++-------------------------- test/toeplitzhankeltests.jl | 33 ++++++ test/toeplitzplanstests.jl | 91 +++++++++++--- 6 files changed, 215 insertions(+), 280 deletions(-) diff --git a/Project.toml b/Project.toml index 3de40b67..7b29f896 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.15.11" +version = "0.15.12" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/src/FastTransforms.jl b/src/FastTransforms.jl index f9d99e74..e9bf6fac 100644 --- a/src/FastTransforms.jl +++ b/src/FastTransforms.jl @@ -8,7 +8,7 @@ using FastGaussQuadrature, FillArrays, LinearAlgebra, @reexport using GenericFFT import Base: convert, unsafe_convert, eltype, ndims, adjoint, transpose, show, - *, \, inv, length, size, view, getindex + *, \, inv, length, size, view, getindex, tail, OneTo import Base.GMP: Limb diff --git a/src/toeplitzhankel.jl b/src/toeplitzhankel.jl index 4e0ab3db..fef934f6 100644 --- a/src/toeplitzhankel.jl +++ b/src/toeplitzhankel.jl @@ -16,75 +16,42 @@ so that `L[:,k] = DL*C[:,k]` and `R[:,k] = DR*C[:,k]`. This allows a Cholesky decomposition in 𝒪(K²N) operations and 𝒪(KN) storage, K = log N log ɛ⁻¹. The tuple storage allows plans applied to each dimension. """ -struct ToeplitzHankelPlan{S, N, M, N1, TP<:ToeplitzPlan{S,N1}} <: Plan{S} - T::NTuple{M,TP} - L::NTuple{M,Matrix{S}} - R::NTuple{M,Matrix{S}} - tmp::Array{S,N1} - dims::NTuple{M,Int} - function ToeplitzHankelPlan{S,N,M,N1,TP}(T::NTuple{M,TP}, L, R, dims) where {S,TP,N,N1,M} +struct ToeplitzHankelPlan{S, N, N1, LowR, TP, Dims} <: Plan{S} + T::TP # A length M Vector or Tuple of ToeplitzPlan + L::LowR # A length M Vector or Tuple of Matrices storing low rank factors of L + R::LowR # A length M Vector or Tuple of Matrices storing low rank factors of R + tmp::Array{S,N1} # A larger dimensional array to transform each scaled array all-at-once + dims::Dims # A length M Vector or Tuple of Int storing the dimensions acted on + function ToeplitzHankelPlan{S,N,N1,LowR,TP,Dims}(T::TP, L::LowR, R::LowR, dims) where {S,N,N1,LowR,TP,Dims} tmp = Array{S}(undef, max.(size.(T)...)...) - new{S,N,M,N1,TP}(T, L, R, tmp, dims) + new{S,N,N1,LowR,TP,Dims}(T, L, R, tmp, dims) end - ToeplitzHankelPlan{S,N,M,N1,TP}(T::NTuple{M,TP}, L, R, dims::Int) where {S,TP,N,N1,M} = - ToeplitzHankelPlan{S,N,M,N1,TP}(T, L, R, (dims,)) end -ToeplitzHankelPlan(T::ToeplitzPlan{S,2}, L::Matrix, R::Matrix, dims=1) where S = - ToeplitzHankelPlan{S, 1, 1, 2, typeof(T)}((T,), (L,), (R,), dims) -ToeplitzHankelPlan(T::ToeplitzPlan{S,3}, L::Matrix, R::Matrix, dims) where S = - ToeplitzHankelPlan{S, 2, 1,3, typeof(T)}((T,), (L,), (R,), dims) +ToeplitzHankelPlan{S,N,M}(T::TP, L::LowR, R::LowR, dims::Dims) where {S,N,M,LowR,TP,Dims} = ToeplitzHankelPlan{S,N,M,LowR,TP,Dims}(T, L, R, dims) +ToeplitzHankelPlan{S,N}(T, L, R, dims) where {S,N} = ToeplitzHankelPlan{S,N,N+1}(T, L, R, dims) +ToeplitzHankelPlan(T::ToeplitzPlan{S,M}, L::Matrix, R::Matrix, dims=1) where {S,M} = ToeplitzHankelPlan{S,M-1,M}((T,), (L,), (R,), dims) -ToeplitzHankelPlan(T::NTuple{2,TP}, L::Tuple, R::Tuple, dims) where {S,TP<:ToeplitzPlan{S,3}} = - ToeplitzHankelPlan{S, 2,2,3, TP}(T, L, R, dims) - -function *(P::ToeplitzHankelPlan{<:Any,1}, v::AbstractVector) - (R,),(L,),(T,),tmp = P.R,P.L,P.T,P.tmp - tmp .= R .* v - T * tmp - tmp .= L .* tmp - sum!(v, tmp) -end - -function _th_applymul1!(v, T, L, R, tmp) - N = size(R,2) - m,n = size(v) - tmp[1:m,1:n,1:N] .= reshape(R,size(R,1),1,N) .* v - T * view(tmp,1:m,1:n,1:N) - view(tmp,1:m,1:n,1:N) .*= reshape(L,size(L,1),1,N) - sum!(v, view(tmp,1:m,1:n,1:N)) -end - -function _th_applymul2!(v, T, L, R, tmp) - N = size(R,2) - m,n = size(v) - tmp[1:m,1:n,1:N] .= reshape(R,1,size(R,1),N) .* v - T * view(tmp,1:m,1:n,1:N) - view(tmp,1:m,1:n,1:N) .*= reshape(L,1,size(L,1),N) - sum!(v, view(tmp,1:m,1:n,1:N)) +_reshape_broadcast(d, R, ::Val{N}, M) where N = reshape(R,ntuple(k -> k == d ? size(R,1) : 1, Val(N))...,M) +function _th_applymul!(d, v::AbstractArray{<:Any,N}, T, L, R, tmp) where N + M = size(R,2) + ax = (axes(v)..., OneTo(M)) + tmp[ax...] .= _reshape_broadcast(d, R, Val(N), M) .* v + T * view(tmp, ax...) + view(tmp,ax...) .*= _reshape_broadcast(d, L, Val(N), M) + sum!(v, view(tmp,ax...)) end -function *(P::ToeplitzHankelPlan{<:Any,2,1}, v::AbstractMatrix) - (R,),(L,),(T,),tmp = P.R,P.L,P.T,P.tmp - if P.dims == (1,) - _th_applymul1!(v, T, L, R, tmp) - else - _th_applymul2!(v, T, L, R, tmp) +function *(P::ToeplitzHankelPlan{<:Any,N}, v::AbstractArray{<:Any,N}) where N + for (R,L,T,d) in zip(P.R,P.L,P.T,P.dims) + _th_applymul!(d, v, T, L, R, P.tmp) end v end -function *(P::ToeplitzHankelPlan{<:Any,2,2}, v::AbstractMatrix) - (R1,R2),(L1,L2),(T1,T2),tmp = P.R,P.L,P.T,P.tmp - - _th_applymul1!(v, T1, L1, R1, tmp) - _th_applymul2!(v, T2, L2, R2, tmp) - - v -end # partial cholesky for a Hankel matrix @@ -166,9 +133,9 @@ function *(P::ChebyshevToLegendrePlanTH, v::AbstractVector{S}) where S v end -function _cheb2leg_rescale1!(V::AbstractMatrix{S}) where S - m,n = size(V) - for j = 1:n +function _cheb2leg_rescale1!(V::AbstractArray{S}) where S + m = size(V,1) + for j = CartesianIndices(tail(axes(V))) ret = zero(S) @inbounds for k = 1:2:m ret += -V[k,j]/(k*(k-2)) @@ -178,24 +145,15 @@ function _cheb2leg_rescale1!(V::AbstractMatrix{S}) where S V end +_dropfirstdim(d::Int) = () +_dropfirstdim(d::Int, m, szs...) = ((d == 1 ? 2 : 1):m, _dropfirstdim(d-1, szs...)...) -function *(P::ChebyshevToLegendrePlanTH, V::AbstractMatrix) +function *(P::ChebyshevToLegendrePlanTH, V::AbstractArray{<:Any,N}) where N m,n = size(V) - dims = P.toeplitzhankel.dims - if dims == (1,) - _cheb2leg_rescale1!(V) - P.toeplitzhankel*view(V,2:m,:) - elseif dims == (2,) - _cheb2leg_rescale1!(transpose(V)) - P.toeplitzhankel*view(V,:,2:n) - else - @assert dims == (1,2) - (R1,R2),(L1,L2),(T1,T2),tmp = P.toeplitzhankel.R,P.toeplitzhankel.L,P.toeplitzhankel.T,P.toeplitzhankel.tmp - - _cheb2leg_rescale1!(V) - _th_applymul1!(view(V,2:m,:), T1, L1, R1, tmp) - _cheb2leg_rescale1!(transpose(V)) - _th_applymul2!(view(V,:,2:n), T2, L2, R2, tmp) + tmp = P.toeplitzhankel.tmp + for (d,R,L,T) in zip(P.toeplitzhankel.dims,P.toeplitzhankel.R,P.toeplitzhankel.L,P.toeplitzhankel.T) + _cheb2leg_rescale1!(PermutedDimsArray(V, _permfirst(d, N))) + _th_applymul!(d, view(V, _dropfirstdim(d, size(V)...)...), T, L, R, tmp) end V end @@ -226,18 +184,14 @@ function _leg2chebuTH_TLC(::Type{S}, mn, d) where {S} (T, (1:n) .* C, C) end - for f in (:leg2cheb, :leg2chebu) plan = Symbol("plan_th_", f, "!") TLC = Symbol("_", f, "TH_TLC") @eval begin - $plan(::Type{S}, mn::Tuple, dims::Int) where {S} = ToeplitzHankelPlan($TLC(S, mn, dims)..., dims) - - function $plan(::Type{S}, mn::NTuple{2,Int}, dims::NTuple{2,Int}) where {S} - @assert dims == (1,2) - T1,L1,C1 = $TLC(S, mn, 1) - T2,L2,C2 = $TLC(S, mn, 2) - ToeplitzHankelPlan((T1,T2), (L1,L2), (C1,C2), dims) + $plan(::Type{S}, mn::NTuple{N,Int}, dims::Int) where {S,N} = ToeplitzHankelPlan($TLC(S, mn, dims)..., dims) + function $plan(::Type{S}, mn::NTuple{N,Int}, dims) where {S,N} + TLCs = $TLC.(S, Ref(mn), dims) + ToeplitzHankelPlan{S,N}(map(first, TLCs), map(TLC -> TLC[2], TLCs), map(last, TLCs), dims) end end end @@ -265,13 +219,11 @@ function _cheb2legTH_TLC(::Type{S}, mn, d) where S T, DL .* C, DR .* C end -plan_th_cheb2leg!(::Type{S}, mn::Tuple, dims::Int) where {S} = ChebyshevToLegendrePlanTH(ToeplitzHankelPlan(_cheb2legTH_TLC(S, mn, dims)..., dims)) +plan_th_cheb2leg!(::Type{S}, mn::NTuple{N,Int}, dims::Int) where {S,N} = ChebyshevToLegendrePlanTH(ToeplitzHankelPlan(_cheb2legTH_TLC(S, mn, dims)..., dims)) -function plan_th_cheb2leg!(::Type{S}, mn::NTuple{2,Int}, dims::NTuple{2,Int}) where {S} - @assert dims == (1,2) - T1,L1,C1 = _cheb2legTH_TLC(S, mn, 1) - T2,L2,C2 = _cheb2legTH_TLC(S, mn, 2) - ChebyshevToLegendrePlanTH(ToeplitzHankelPlan((T1,T2), (L1,L2), (C1,C2), dims)) +function plan_th_cheb2leg!(::Type{S}, mn::NTuple{N,Int}, dims) where {S,N} + TLCs = _cheb2legTH_TLC.(S, Ref(mn), dims) + ChebyshevToLegendrePlanTH(ToeplitzHankelPlan{S,N}(map(first, TLCs), map(TLC -> TLC[2], TLCs), map(last, TLCs), dims)) end @@ -337,7 +289,7 @@ _good_plan_th_ultra2ultra!(::Type{S}, mn, λ₁, λ₂, dims::Int) where S = Toe function _good_plan_th_ultra2ultra!(::Type{S}, mn::NTuple{2,Int}, λ₁, λ₂, dims::NTuple{2,Int}) where S T1,L1,C1 = _ultra2ultraTH_TLC(S, mn, λ₁, λ₂, 1) T2,L2,C2 = _ultra2ultraTH_TLC(S, mn, λ₁, λ₂, 2) - ToeplitzHankelPlan((T1,T2), (L1,L2), (C1,C2), dims) + ToeplitzHankelPlan{S,2}((T1,T2), (L1,L2), (C1,C2), dims) end @@ -515,7 +467,7 @@ _good_plan_th_jac2jac!(::Type{S}, mn, α, β, γ, δ, dims::Int) where S = Toepl function _good_plan_th_jac2jac!(::Type{S}, mn::NTuple{2,Int}, α, β, γ, δ, dims::NTuple{2,Int}) where S T1,L1,C1 = _jac2jacTH_TLC(S, mn, α, β, γ, δ, 1) T2,L2,C2 = _jac2jacTH_TLC(S, mn, α, β, γ, δ, 2) - ToeplitzHankelPlan((T1,T2), (L1,L2), (C1,C2), dims) + ToeplitzHankelPlan{S,2}((T1,T2), (L1,L2), (C1,C2), dims) end @@ -685,10 +637,8 @@ end for f in (:th_leg2cheb, :th_cheb2leg, :th_leg2chebu) plan = Symbol("plan_", f, "!") @eval begin - $plan(::Type{S}, mn::NTuple{N,Int}, dims::UnitRange) where {N,S} = $plan(S, mn, tuple(dims...)) - $plan(::Type{S}, mn::Tuple{Int}, dims::Tuple{Int}=(1,)) where {S} = $plan(S, mn, dims...) - $plan(::Type{S}, (m,n)::NTuple{2,Int}) where {S} = $plan(S, (m,n), (1,2)) $plan(arr::AbstractArray{T}, dims...) where T = $plan(T, size(arr), dims...) + $plan(::Type{S}, mn::NTuple{N,Int}) where {S,N} = $plan(S, mn, ntuple(identity,Val(N))) $f(v, dims...) = $plan(eltype(v), size(v), dims...)*copy(v) end end diff --git a/src/toeplitzplans.jl b/src/toeplitzplans.jl index ab08d1f7..42d24062 100644 --- a/src/toeplitzplans.jl +++ b/src/toeplitzplans.jl @@ -1,39 +1,46 @@ using FFTW import FFTW: plan_r2r! -struct ToeplitzPlan{T, N, M, S, VECS<:Tuple{Vararg{Vector{S}}}, P<:Plan{S}, Pi<:Plan{S}} <: Plan{T} - vectors::VECS + +""" + ToeplitzPlan + +applies Toeplitz matrices fast along each dimension. +""" + +struct ToeplitzPlan{T, N, Dims, S, VECS, P<:Plan{S}, Pi<:Plan{S}} <: Plan{T} + vectors::VECS # Vector or Tuple of storage tmp::Array{S,N} dft::P idft::Pi - dims::NTuple{M,Int} + dims::Dims end -ToeplitzPlan{T}(v::AbstractVector, tmp, dft, idft, dims) where T = ToeplitzPlan{T}((v,), tmp, dft, idft, dims) -ToeplitzPlan{T}(v::Tuple{Vararg{Vector{S}}}, tmp::Array{S,N}, dft::Plan{S}, idft::Plan{S}, dims::NTuple{M,Int}) where {T,S,N,M} = ToeplitzPlan{T,N,M,S,typeof(v),typeof(dft), typeof(idft)}(v, tmp, dft, idft, dims) -ToeplitzPlan{T}(v::Tuple{Vararg{Vector{S}}}, tmp::Array{S,N}, dft::Plan{S}, idft::Plan{S}, dims::Int) where {T,S,N} = ToeplitzPlan{T}(v, tmp, dft, idft, (dims,)) +ToeplitzPlan{T}(v, tmp::Array{S,N}, dft::Plan{S}, idft::Plan{S}, dims) where {T,S,N} = ToeplitzPlan{T,N,typeof(dims),S,typeof(v),typeof(dft), typeof(idft)}(v, tmp, dft, idft, dims) + -size(A::ToeplitzPlan{<:Any,1}) = ((length(A.tmp)+1) ÷ 2,) -function size(A::ToeplitzPlan{<:Any,2,1}) - if A.dims == (1,) - ((size(A.tmp,1)+1) ÷ 2, size(A.tmp,2)) - else # A.dims == (2,) - (size(A.tmp,1), (size(A.tmp,2)+1) ÷ 2) +divdimby2(d::Int, sz1, szs...) = isone(d) ? ((sz1 + 1) ÷ 2, szs...) : (sz1, divdimby2(d-1, szs...)...) +muldimby2(d::Int, sz1, szs...) = isone(d) ? (max(0,2sz1 - 1), szs...) : (sz1, muldimby2(d-1, szs...)...) + +function toeplitzplan_size(dims, szs) + ret = szs + for d in dims + ret = divdimby2(d, ret...) end + ret end -function size(A::ToeplitzPlan{<:Any,3,1}) - if A.dims == (1,) - ((size(A.tmp,1)+1) ÷ 2, size(A.tmp,2), size(A.tmp,3)) - elseif A.dims == (2,) - (size(A.tmp,1), (size(A.tmp,2)+1) ÷ 2, size(A.tmp,3)) - else - (size(A.tmp,1), size(A.tmp,2), (size(A.tmp,3)+1) ÷ 2) +function to_toeplitzplan_size(dims, szs) + ret = szs + for d in dims + ret = muldimby2(d, ret...) end + ret end -size(A::ToeplitzPlan{<:Any,2,2}) = ((size(A.tmp,1)+1) ÷ 2, (size(A.tmp,2)+1) ÷ 2) +size(A::ToeplitzPlan) = toeplitzplan_size(A.dims, size(A.tmp)) + # based on ToeplitzMatrices.jl """ @@ -44,101 +51,31 @@ Return real-valued part of `x` if `T` is a type of a real number, and `x` otherw maybereal(::Type, x) = x maybereal(::Type{<:Real}, x) = real(x) -function *(A::ToeplitzPlan{T,1}, x::AbstractVector{T}) where T - vc,tmp,dft,idft = A.vectors[1],A.tmp, A.dft,A.idft - S = eltype(tmp) - N = length(tmp) - n = length(x) - if 2n-1 ≠ N - throw(DimensionMismatch("Toeplitz plan does not match size of input")) - end - copyto!(view(tmp, 1:n), x) - fill!(view(tmp, n+1:N), zero(S)) - dft * tmp - tmp .*= vc - idft * tmp - @inbounds for k = 1:n - x[k] = maybereal(T, tmp[k]) - end - x -end +function *(A::ToeplitzPlan{T,N}, X::AbstractArray{T,N}) where {T,N} + vcs,Y,dft,idft,dims = A.vectors,A.tmp, A.dft,A.idft,A.dims -function *(A::ToeplitzPlan{T,2,1, S}, x::AbstractMatrix{T}) where {T,S} - vc,tmp,dft,idft = A.vectors[1],A.tmp, A.dft, A.idft - M,N = size(tmp) - m,n = size(x) + isempty(X) && return X - if isempty(x) - return x - end + fill!(Y, zero(eltype(Y))) + copyto!(view(Y, axes(X)...), X) - if A.dims == (1,) - copyto!(view(tmp, 1:m, :), x) - fill!(view(tmp, m+1:M, :), zero(S)) - if !isempty(tmp) - dft * tmp - end - tmp .= vc .* tmp - else - @assert A.dims == (2,) - copyto!(view(tmp, :, 1:n), x) - fill!(view(tmp, :, n+1:N), zero(S)) - dft * tmp - tmp .= tmp .* transpose(vc) + # Fourier transform each dimension + dft * Y + + # Multiply by a diagonal matrix along each dimension by permuting + # to first dimension + for (vc,d) in zip(vcs,dims) + Ỹ = PermutedDimsArray(Y, _permfirst(d, N)) + Ỹ .= vc .* Ỹ end - idft * tmp - x .= maybereal.(T, view(tmp,1:m,1:n)) -end + # Transform back + idft * Y -function *(A::ToeplitzPlan{T,2,2, S}, X::AbstractMatrix{T}) where {T,S} - vcs,tmp,dft,idft = A.vectors,A.tmp, A.dft,A.idft - vc1,vc2 = vcs - M,N = size(tmp) - m,n = size(X) - - @assert A.dims == (1,2) - copyto!(view(tmp, 1:m, 1:n), X) - fill!(view(tmp, m+1:M, :), zero(S)) - fill!(view(tmp, 1:m, n+1:N), zero(S)) - dft * tmp - tmp .= vc1 .* tmp .* transpose(vc2) - idft * tmp - @inbounds for k = 1:m, j = 1:n - X[k,j] = maybereal(T, tmp[k,j]) - end + X .= maybereal.(T, view(Y, axes(X)...)) X end -function *(A::ToeplitzPlan{T,3,1, S}, x::AbstractArray{T,3}) where {T,S} - vc,tmp,dft,idft = A.vectors[1],A.tmp, A.dft,A.idft - M,N,L = size(tmp) - m,n,l = size(x) - - if A.dims == (1,) - copyto!(view(tmp, 1:m, :, :), x) - fill!(view(tmp, m+1:M, :, :), zero(S)) - dft * tmp - tmp .= vc .* tmp - elseif A.dims == (2,) - copyto!(view(tmp, :, 1:n, :), x) - fill!(view(tmp, :, n+1:N, :), zero(S)) - dft * tmp - tmp .= tmp .* transpose(vc) - else - copyto!(view(tmp, :, :, 1:l), x) - fill!(view(tmp, :, :, l+1:L), zero(S)) - dft * tmp - tmp .= tmp .* reshape(vc, 1, 1, L) - end - idft * tmp - @inbounds for k = 1:m, j = 1:n, ℓ = 1:l - x[k,j,ℓ] = maybereal(T, tmp[k,j,ℓ]) - end - x -end - - function uppertoeplitz_padvec(v::AbstractVector{T}) where T n = length(v) @@ -151,73 +88,25 @@ function uppertoeplitz_padvec(v::AbstractVector{T}) where T tmp end -function plan_uppertoeplitz!(v::AbstractVector{T}) where T - tmp = uppertoeplitz_padvec(v) - dft = plan_fft!(tmp) - idft = plan_ifft!(similar(tmp)) - return ToeplitzPlan{float(T)}(dft * tmp, similar(tmp), dft, idft, (1,)) -end +safe_fft!(A) = isempty(A) ? A : fft!(A) -# TODO: support different transforms -# function plan_uppertoeplitz!(v1::AbstractVector{T}, v2::AbstractVector{T}) where T -# S = float(T) -# m,n = length(v1), length(v2) -# tmp = zeros(S, 2m-1, 2n-1) -# pv1 = uppertoeplitz_padvec(v1) -# pv2 = uppertoeplitz_padvec(v2) -# dft = plan_r2r!(tmp, FFTW.R2HC) -# return ToeplitzPlan((r2r!(pv1, FFTW.R2HC), r2r!(pv2, FFTW.R2HC)), tmp, dft, 1:2) -# end - -function plan_uppertoeplitz!(v::AbstractVector{T}, szs::NTuple{2,Int}, dim::Int) where T - S = complex(float(T)) - m,n = szs - if isone(dim) - tmp = zeros(S, max(0,2m-1), n) - pv = uppertoeplitz_padvec(v[1:m]) - else # dim == 2 - tmp = zeros(S, m, max(0,2n-1)) - pv = uppertoeplitz_padvec(v[1:n]) - end - if isempty(tmp) - # dummy plans just to create type - dft = plan_fft!(similar(tmp, 1, 1), dim) - idft = plan_ifft!(similar(tmp, 1, 1), dim) - ToeplitzPlan{float(T)}(pv, tmp, dft, idft, dim) - else - dft = plan_fft!(tmp, dim) - idft = plan_ifft!(similar(tmp), dim) - return ToeplitzPlan{float(T)}(fft!(pv), tmp, dft, idft, dim) - end -end +uppertoeplitz_vecs(v, dims::AbstractVector, szs) = [safe_fft!(uppertoeplitz_padvec(v[1:szs[d]])) for d in dims] +uppertoeplitz_vecs(v, dims::Tuple{}, szs) = () +uppertoeplitz_vecs(v, dims::Tuple, szs) = (safe_fft!(uppertoeplitz_padvec(v[1:szs[first(dims)]])), uppertoeplitz_vecs(v, tail(dims), szs)...) +uppertoeplitz_vecs(v, d::Int, szs) = (safe_fft!(uppertoeplitz_padvec(v[1:szs[d]])),) -function plan_uppertoeplitz!(v::AbstractVector{T}, szs::NTuple{3,Int}, dim::Int) where T - S = complex(float(T)) - m,n,l = szs - if isone(dim) - tmp = zeros(S, 2m-1, n, l) - pv = uppertoeplitz_padvec(v[1:m]) - elseif dim == 2 - tmp = zeros(S, m, 2n-1, l) - pv = uppertoeplitz_padvec(v[1:n]) - else - @assert dim == 3 - tmp = zeros(S, m, n, 2l-1) - pv = uppertoeplitz_padvec(v[1:l]) - end - dft = plan_fft!(tmp, dim) - idft = plan_ifft!(similar(tmp), dim) - return ToeplitzPlan{float(T)}(fft!(pv), tmp, dft, idft, dim) -end -function plan_uppertoeplitz!(v::AbstractVector{T}, szs::NTuple{2,Int}, dim=(1,2)) where T - @assert dim == (1,2) +# allow FFT to work by making sure tmp is non-empty +safe_tmp(tmp::AbstractArray{<:Any,N}) where N = isempty(tmp) ? similar(tmp, ntuple(_ -> 1, Val(N))...) : tmp + +function plan_uppertoeplitz!(v::AbstractVector{T}, szs::NTuple{N,Int}, dim=ntuple(identity,Val(N))) where {T,N} S = complex(float(T)) - m,n = szs - tmp = zeros(S, 2m-1, 2n-1) - pv1 = uppertoeplitz_padvec(v[1:m]) - pv2 = uppertoeplitz_padvec(v[1:n]) - dft = plan_fft!(tmp, dim) - idft = plan_ifft!(similar(tmp), dim) - return ToeplitzPlan{float(T)}((fft!(pv1), fft!(pv2)), tmp, dft, idft, dim) + + tmp = zeros(S, to_toeplitzplan_size(dim, szs)...) + dft = plan_fft!(safe_tmp(tmp), dim) + idft = plan_ifft!(safe_tmp(similar(tmp)), dim) + + return ToeplitzPlan{float(T)}(uppertoeplitz_vecs(v, dim, szs), tmp, dft, idft, dim) end + +plan_uppertoeplitz!(v::AbstractVector{T}) where T = plan_uppertoeplitz!(v, size(v)) diff --git a/test/toeplitzhankeltests.jl b/test/toeplitzhankeltests.jl index 19274747..b90ce94b 100644 --- a/test/toeplitzhankeltests.jl +++ b/test/toeplitzhankeltests.jl @@ -130,4 +130,37 @@ Random.seed!(0) @test norm(v - th_cheb2leg(th_leg2cheb(v)), Inf) ≤ 1E-13 @test norm(v - th_cheb2leg(th_leg2cheb(v)))/norm(v) ≤ 1E-14 end + + @testset "tensor" begin + X = randn(5,4,3) + for trans in (th_leg2cheb, th_cheb2leg) + Y = trans(X, 1) + for ℓ = 1:size(X,3) + @test Y[:,:,ℓ] ≈ trans(X[:,:,ℓ],1) + end + Y = trans(X, 2) + for ℓ = 1:size(X,3) + @test Y[:,:,ℓ] ≈ trans(X[:,:,ℓ],2) + end + Y = trans(X, 3) + for j = 1:size(X,2) + @test Y[:,j,:] ≈ trans(X[:,j,:],2) + end + + Y = trans(X, (1,3)) + for j = 1:size(X,2) + @test Y[:,j,:] ≈ trans(X[:,j,:]) + end + + Y = trans(X, 1:3) + M = copy(X) + for j = 1:size(X,3) + M[:,:,j] = trans(M[:,:,j]) + end + for k = 1:size(X,1), j=1:size(X,2) + M[k,j,:] = trans(M[k,j,:]) + end + @test M ≈ Y + end + end end \ No newline at end of file diff --git a/test/toeplitzplanstests.jl b/test/toeplitzplanstests.jl index e56d8c3e..6ea6a095 100644 --- a/test/toeplitzplanstests.jl +++ b/test/toeplitzplanstests.jl @@ -34,23 +34,86 @@ import FastTransforms: plan_uppertoeplitz! @testset "Tensor" begin T = [1 2 3; 0 1 2; 0 0 1] - X = randn(3,3,3) - P = plan_uppertoeplitz!([1,2,3], size(X), 1) - PX = P * copy(X) - for ℓ = 1:size(X,3) - @test PX[:,:,ℓ] ≈ T*X[:,:,ℓ] - end + @testset "3D" begin + X = randn(3,3,3) + P = plan_uppertoeplitz!([1,2,3], size(X), 1) + PX = P * copy(X) + for ℓ = 1:size(X,3) + @test PX[:,:,ℓ] ≈ T*X[:,:,ℓ] + end - P = plan_uppertoeplitz!([1,2,3], size(X), 2) - PX = P * copy(X) - for ℓ = 1:size(X,3) - @test PX[:,:,ℓ] ≈ X[:,:,ℓ]*T' + P = plan_uppertoeplitz!([1,2,3], size(X), 2) + PX = P * copy(X) + for ℓ = 1:size(X,3) + @test PX[:,:,ℓ] ≈ X[:,:,ℓ]*T' + end + + P = plan_uppertoeplitz!([1,2,3], size(X), 3) + PX = P * copy(X) + for j = 1:size(X,2) + @test PX[:,j,:] ≈ X[:,j,:]*T' + end + + P = plan_uppertoeplitz!([1,2,3], size(X), (1,3)) + PX = P * copy(X) + for j = 1:size(X,2) + @test PX[:,j,:] ≈ T*X[:,j,:]*T' + end + + P = plan_uppertoeplitz!([1,2,3], size(X), 1:3) + PX = P * copy(X) + M = copy(X) + for j = 1:size(X,3) + M[:,:,j] = T*M[:,:,j]*T' + end + for k = 1:size(X,1) + M[k,:,:] = M[k,:,:]*T' + end + @test M ≈ PX end - P = plan_uppertoeplitz!([1,2,3], size(X), 3) - PX = P * copy(X) - for j = 1:size(X,2) - @test PX[:,j,:] ≈ X[:,j,:]*T' + @testset "4D" begin + X = randn(3,3,3,3) + P = plan_uppertoeplitz!([1,2,3], size(X), 1) + PX = P * copy(X) + for ℓ = 1:size(X,3), m = 1:size(X,4) + @test PX[:,:,ℓ,m] ≈ T*X[:,:,ℓ,m] + end + + P = plan_uppertoeplitz!([1,2,3], size(X), 2) + PX = P * copy(X) + for ℓ = 1:size(X,3), m = 1:size(X,4) + @test PX[:,:,ℓ,m] ≈ X[:,:,ℓ,m]*T' + end + + P = plan_uppertoeplitz!([1,2,3], size(X), 3) + PX = P * copy(X) + for j = 1:size(X,2), m = 1:size(X,4) + @test PX[:,j,:,m] ≈ X[:,j,:,m]*T' + end + + P = plan_uppertoeplitz!([1,2,3], size(X), 4) + PX = P * copy(X) + for k = 1:size(X,1), j = 1:size(X,2) + @test PX[k,j,:,:] ≈ X[k,j,:,:]*T' + end + + P = plan_uppertoeplitz!([1,2,3], size(X), (1,3)) + PX = P * copy(X) + for j = 1:size(X,2), m=1:size(X,4) + @test PX[:,j,:,m] ≈ T*X[:,j,:,m]*T' + end + + P = plan_uppertoeplitz!([1,2,3], size(X), 1:4) + PX = P * copy(X) + M = copy(X) + for ℓ = 1:size(X,3), m = 1:size(X,4) + M[:,:,ℓ,m] = T*M[:,:,ℓ,m]*T' + end + for k = 1:size(X,1), j = 1:size(X,2) + M[k,j,:,:] = T*M[k,j,:,:]*T' + end + @test M ≈ PX end end From 90d1daefccff78ac9f45ee0d1ba487d0f1915569 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 13 Nov 2023 14:41:38 +0000 Subject: [PATCH 10/25] Cheb2LegTH <: Plan (#234) * Cheb2LegTH <: Plan * avoid stack overflow in different sizes --- src/chebyshevtransform.jl | 4 ++++ src/toeplitzhankel.jl | 2 +- test/chebyshevtests.jl | 5 +++++ 3 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src/chebyshevtransform.jl b/src/chebyshevtransform.jl index 01dd61af..541d9933 100644 --- a/src/chebyshevtransform.jl +++ b/src/chebyshevtransform.jl @@ -2,6 +2,8 @@ abstract type ChebyshevPlan{T} <: Plan{T} end +*(P::ChebyshevPlan{T}, x::AbstractArray{T}) where T = error("Plan applied to wrong size array") + size(P::ChebyshevPlan) = isdefined(P, :plan) ? size(P.plan) : (0,) length(P::ChebyshevPlan) = isdefined(P, :plan) ? length(P.plan) : 0 @@ -151,6 +153,8 @@ end ldiv!(_prod_size(size(y), d), y) end + + function *(P::ChebyshevTransformPlan{T,1,K,true,N}, x::AbstractArray{T,N}) where {T,K,N} isempty(x) && return x diff --git a/src/toeplitzhankel.jl b/src/toeplitzhankel.jl index fef934f6..383ff9ae 100644 --- a/src/toeplitzhankel.jl +++ b/src/toeplitzhankel.jl @@ -118,7 +118,7 @@ end -struct ChebyshevToLegendrePlanTH{TH} +struct ChebyshevToLegendrePlanTH{S,TH<:ToeplitzHankelPlan{S}} <: Plan{S} toeplitzhankel::TH end diff --git a/test/chebyshevtests.jl b/test/chebyshevtests.jl index d6f704bf..051d5f51 100644 --- a/test/chebyshevtests.jl +++ b/test/chebyshevtests.jl @@ -462,4 +462,9 @@ using FastTransforms, Test @test F \ (F*X) ≈ F * (F\X) ≈ X end end + + @testset "incompatible shapes" begin + @test_throws ErrorException plan_chebyshevtransform(randn(5)) * randn(5,5) + @test_throws ErrorException plan_ichebyshevtransform(randn(5)) * randn(5,5) + end end From e1db1ca2c4f3516ff3ac7beaa80cec197a84444d Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Tue, 14 Nov 2023 14:07:23 +0000 Subject: [PATCH 11/25] 4D Chebyshev transform (#235) --- Project.toml | 2 +- src/chebyshevtransform.jl | 84 ++++------------- test/chebyshevtests.jl | 183 ++++++++++++++++++++------------------ 3 files changed, 113 insertions(+), 156 deletions(-) diff --git a/Project.toml b/Project.toml index 7b29f896..3d45e5c3 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.15.12" +version = "0.15.13" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/src/chebyshevtransform.jl b/src/chebyshevtransform.jl index 541d9933..d55b52c3 100644 --- a/src/chebyshevtransform.jl +++ b/src/chebyshevtransform.jl @@ -58,76 +58,22 @@ plan_chebyshevtransform(x::AbstractArray, dims...; kws...) = plan_chebyshevtrans @inline _plan_mul!(y::AbstractArray{T}, P::Plan{T}, x::AbstractArray) where T = mul!(y, P, convert(Array{T}, x)) +for op in (:ldiv, :lmul) + op_dim_begin! = Symbol(string(op) * "_dim_begin!") + op_dim_end! = Symbol(string(op) * "_dim_end!") + op! = Symbol(string(op) * "!") + @eval begin + function $op_dim_begin!(α, d::Number, y::AbstractArray{<:Any,N}) where N + # scale just the d-th dimension by permuting it to the first + ỹ = PermutedDimsArray(y, _permfirst(d, N)) + $op!(α, view(ỹ, 1, ntuple(_ -> :, Val(N-1))...)) + end -ldiv_dim_begin!(α, d::Number, y::AbstractVector) = y[1] /= α -function ldiv_dim_begin!(α, d::Number, y::AbstractMatrix) - if isone(d) - ldiv!(α, @view(y[1,:])) - else - ldiv!(α, @view(y[:,1])) - end -end -function ldiv_dim_begin!(α, d::Number, y::AbstractArray{<:Any,3}) - if isone(d) - ldiv!(α, @view(y[1,:,:])) - elseif d == 2 - ldiv!(α, @view(y[:,1,:])) - else # d == 3 - ldiv!(α, @view(y[:,:,1])) - end -end - -ldiv_dim_end!(α, d::Number, y::AbstractVector) = y[end] /= α -function ldiv_dim_end!(α, d::Number, y::AbstractMatrix) - if isone(d) - ldiv!(α, @view(y[end,:])) - else - ldiv!(α, @view(y[:,end])) - end -end -function ldiv_dim_end!(α, d::Number, y::AbstractArray{<:Any,3}) - if isone(d) - ldiv!(α, @view(y[end,:,:])) - elseif d == 2 - ldiv!(α, @view(y[:,end,:])) - else # d == 3 - ldiv!(α, @view(y[:,:,end])) - end -end - -lmul_dim_begin!(α, d::Number, y::AbstractVector) = y[1] *= α -function lmul_dim_begin!(α, d::Number, y::AbstractMatrix) - if isone(d) - lmul!(α, @view(y[1,:])) - else - lmul!(α, @view(y[:,1])) - end -end -function lmul_dim_begin!(α, d::Number, y::AbstractArray{<:Any,3}) - if isone(d) - lmul!(α, @view(y[1,:,:])) - elseif d == 2 - lmul!(α, @view(y[:,1,:])) - else # d == 3 - lmul!(α, @view(y[:,:,1])) - end -end - -lmul_dim_end!(α, d::Number, y::AbstractVector) = y[end] *= α -function lmul_dim_end!(α, d::Number, y::AbstractMatrix) - if isone(d) - lmul!(α, @view(y[end,:])) - else - lmul!(α, @view(y[:,end])) - end -end -function lmul_dim_end!(α, d::Number, y::AbstractArray{<:Any,3}) - if isone(d) - lmul!(α, @view(y[end,:,:])) - elseif d == 2 - lmul!(α, @view(y[:,end,:])) - else # d == 3 - lmul!(α, @view(y[:,:,end])) + function $op_dim_end!(α, d::Number, y::AbstractArray{<:Any,N}) where N + # scale just the d-th dimension by permuting it to the first + ỹ = PermutedDimsArray(y, _permfirst(d, N)) + $op!(α, view(ỹ, size(ỹ,1), ntuple(_ -> :, Val(N-1))...)) + end end end diff --git a/test/chebyshevtests.jl b/test/chebyshevtests.jl index 051d5f51..2a82978c 100644 --- a/test/chebyshevtests.jl +++ b/test/chebyshevtests.jl @@ -322,96 +322,107 @@ using FastTransforms, Test end @testset "tensor" begin - X = randn(4,5,6) - X̃ = similar(X) - @testset "chebyshevtransform" begin - for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = chebyshevtransform(X[:,k,j]) end - @test @inferred(chebyshevtransform(X,1)) ≈ @inferred(chebyshevtransform!(copy(X),1)) ≈ X̃ - for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = chebyshevtransform(X[k,:,j]) end - @test chebyshevtransform(X,2) ≈ chebyshevtransform!(copy(X),2) ≈ X̃ - for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = chebyshevtransform(X[k,j,:]) end - @test chebyshevtransform(X,3) ≈ chebyshevtransform!(copy(X),3) ≈ X̃ - - for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = chebyshevtransform(X[:,k,j],Val(2)) end - @test @inferred(chebyshevtransform(X,Val(2),1)) ≈ @inferred(chebyshevtransform!(copy(X),Val(2),1)) ≈ X̃ - for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = chebyshevtransform(X[k,:,j],Val(2)) end - @test chebyshevtransform(X,Val(2),2) ≈ chebyshevtransform!(copy(X),Val(2),2) ≈ X̃ - for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = chebyshevtransform(X[k,j,:],Val(2)) end - @test chebyshevtransform(X,Val(2),3) ≈ chebyshevtransform!(copy(X),Val(2),3) ≈ X̃ - - @test @inferred(chebyshevtransform(X)) ≈ @inferred(chebyshevtransform!(copy(X))) ≈ chebyshevtransform(chebyshevtransform(chebyshevtransform(X,1),2),3) - @test @inferred(chebyshevtransform(X,Val(2))) ≈ @inferred(chebyshevtransform!(copy(X),Val(2))) ≈ chebyshevtransform(chebyshevtransform(chebyshevtransform(X,Val(2),1),Val(2),2),Val(2),3) - end - - @testset "ichebyshevtransform" begin - for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = ichebyshevtransform(X[:,k,j]) end - @test @inferred(ichebyshevtransform(X,1)) ≈ @inferred(ichebyshevtransform!(copy(X),1)) ≈ X̃ - for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = ichebyshevtransform(X[k,:,j]) end - @test ichebyshevtransform(X,2) ≈ ichebyshevtransform!(copy(X),2) ≈ X̃ - for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = ichebyshevtransform(X[k,j,:]) end - @test ichebyshevtransform(X,3) ≈ ichebyshevtransform!(copy(X),3) ≈ X̃ - - for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = ichebyshevtransform(X[:,k,j],Val(2)) end - @test @inferred(ichebyshevtransform(X,Val(2),1)) ≈ @inferred(ichebyshevtransform!(copy(X),Val(2),1)) ≈ X̃ - for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = ichebyshevtransform(X[k,:,j],Val(2)) end - @test ichebyshevtransform(X,Val(2),2) ≈ ichebyshevtransform!(copy(X),Val(2),2) ≈ X̃ - for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = ichebyshevtransform(X[k,j,:],Val(2)) end - @test ichebyshevtransform(X,Val(2),3) ≈ ichebyshevtransform!(copy(X),Val(2),3) ≈ X̃ - - @test @inferred(ichebyshevtransform(X)) ≈ @inferred(ichebyshevtransform!(copy(X))) ≈ ichebyshevtransform(ichebyshevtransform(ichebyshevtransform(X,1),2),3) - @test @inferred(ichebyshevtransform(X,Val(2))) ≈ @inferred(ichebyshevtransform!(copy(X),Val(2))) ≈ ichebyshevtransform(ichebyshevtransform(ichebyshevtransform(X,Val(2),1),Val(2),2),Val(2),3) - - @test ichebyshevtransform(chebyshevtransform(X)) ≈ X - @test chebyshevtransform(ichebyshevtransform(X)) ≈ X - end - - @testset "chebyshevutransform" begin - for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = chebyshevutransform(X[:,k,j]) end - @test @inferred(chebyshevutransform(X,1)) ≈ @inferred(chebyshevutransform!(copy(X),1)) ≈ X̃ - for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = chebyshevutransform(X[k,:,j]) end - @test chebyshevutransform(X,2) ≈ chebyshevutransform!(copy(X),2) ≈ X̃ - for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = chebyshevutransform(X[k,j,:]) end - @test chebyshevutransform(X,3) ≈ chebyshevutransform!(copy(X),3) ≈ X̃ - - for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = chebyshevutransform(X[:,k,j],Val(2)) end - @test @inferred(chebyshevutransform(X,Val(2),1)) ≈ @inferred(chebyshevutransform!(copy(X),Val(2),1)) ≈ X̃ - for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = chebyshevutransform(X[k,:,j],Val(2)) end - @test chebyshevutransform(X,Val(2),2) ≈ chebyshevutransform!(copy(X),Val(2),2) ≈ X̃ - for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = chebyshevutransform(X[k,j,:],Val(2)) end - @test chebyshevutransform(X,Val(2),3) ≈ chebyshevutransform!(copy(X),Val(2),3) ≈ X̃ - - @test @inferred(chebyshevutransform(X)) ≈ @inferred(chebyshevutransform!(copy(X))) ≈ chebyshevutransform(chebyshevutransform(chebyshevutransform(X,1),2),3) - @test @inferred(chebyshevutransform(X,Val(2))) ≈ @inferred(chebyshevutransform!(copy(X),Val(2))) ≈ chebyshevutransform(chebyshevutransform(chebyshevutransform(X,Val(2),1),Val(2),2),Val(2),3) + @testset "3D" begin + X = randn(4,5,6) + X̃ = similar(X) + @testset "chebyshevtransform" begin + for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = chebyshevtransform(X[:,k,j]) end + @test @inferred(chebyshevtransform(X,1)) ≈ @inferred(chebyshevtransform!(copy(X),1)) ≈ X̃ + for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = chebyshevtransform(X[k,:,j]) end + @test chebyshevtransform(X,2) ≈ chebyshevtransform!(copy(X),2) ≈ X̃ + for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = chebyshevtransform(X[k,j,:]) end + @test chebyshevtransform(X,3) ≈ chebyshevtransform!(copy(X),3) ≈ X̃ + + for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = chebyshevtransform(X[:,k,j],Val(2)) end + @test @inferred(chebyshevtransform(X,Val(2),1)) ≈ @inferred(chebyshevtransform!(copy(X),Val(2),1)) ≈ X̃ + for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = chebyshevtransform(X[k,:,j],Val(2)) end + @test chebyshevtransform(X,Val(2),2) ≈ chebyshevtransform!(copy(X),Val(2),2) ≈ X̃ + for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = chebyshevtransform(X[k,j,:],Val(2)) end + @test chebyshevtransform(X,Val(2),3) ≈ chebyshevtransform!(copy(X),Val(2),3) ≈ X̃ + + @test @inferred(chebyshevtransform(X)) ≈ @inferred(chebyshevtransform!(copy(X))) ≈ chebyshevtransform(chebyshevtransform(chebyshevtransform(X,1),2),3) + @test @inferred(chebyshevtransform(X,Val(2))) ≈ @inferred(chebyshevtransform!(copy(X),Val(2))) ≈ chebyshevtransform(chebyshevtransform(chebyshevtransform(X,Val(2),1),Val(2),2),Val(2),3) + end + + @testset "ichebyshevtransform" begin + for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = ichebyshevtransform(X[:,k,j]) end + @test @inferred(ichebyshevtransform(X,1)) ≈ @inferred(ichebyshevtransform!(copy(X),1)) ≈ X̃ + for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = ichebyshevtransform(X[k,:,j]) end + @test ichebyshevtransform(X,2) ≈ ichebyshevtransform!(copy(X),2) ≈ X̃ + for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = ichebyshevtransform(X[k,j,:]) end + @test ichebyshevtransform(X,3) ≈ ichebyshevtransform!(copy(X),3) ≈ X̃ + + for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = ichebyshevtransform(X[:,k,j],Val(2)) end + @test @inferred(ichebyshevtransform(X,Val(2),1)) ≈ @inferred(ichebyshevtransform!(copy(X),Val(2),1)) ≈ X̃ + for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = ichebyshevtransform(X[k,:,j],Val(2)) end + @test ichebyshevtransform(X,Val(2),2) ≈ ichebyshevtransform!(copy(X),Val(2),2) ≈ X̃ + for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = ichebyshevtransform(X[k,j,:],Val(2)) end + @test ichebyshevtransform(X,Val(2),3) ≈ ichebyshevtransform!(copy(X),Val(2),3) ≈ X̃ + + @test @inferred(ichebyshevtransform(X)) ≈ @inferred(ichebyshevtransform!(copy(X))) ≈ ichebyshevtransform(ichebyshevtransform(ichebyshevtransform(X,1),2),3) + @test @inferred(ichebyshevtransform(X,Val(2))) ≈ @inferred(ichebyshevtransform!(copy(X),Val(2))) ≈ ichebyshevtransform(ichebyshevtransform(ichebyshevtransform(X,Val(2),1),Val(2),2),Val(2),3) + + @test ichebyshevtransform(chebyshevtransform(X)) ≈ X + @test chebyshevtransform(ichebyshevtransform(X)) ≈ X + end + + @testset "chebyshevutransform" begin + for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = chebyshevutransform(X[:,k,j]) end + @test @inferred(chebyshevutransform(X,1)) ≈ @inferred(chebyshevutransform!(copy(X),1)) ≈ X̃ + for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = chebyshevutransform(X[k,:,j]) end + @test chebyshevutransform(X,2) ≈ chebyshevutransform!(copy(X),2) ≈ X̃ + for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = chebyshevutransform(X[k,j,:]) end + @test chebyshevutransform(X,3) ≈ chebyshevutransform!(copy(X),3) ≈ X̃ + + for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = chebyshevutransform(X[:,k,j],Val(2)) end + @test @inferred(chebyshevutransform(X,Val(2),1)) ≈ @inferred(chebyshevutransform!(copy(X),Val(2),1)) ≈ X̃ + for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = chebyshevutransform(X[k,:,j],Val(2)) end + @test chebyshevutransform(X,Val(2),2) ≈ chebyshevutransform!(copy(X),Val(2),2) ≈ X̃ + for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = chebyshevutransform(X[k,j,:],Val(2)) end + @test chebyshevutransform(X,Val(2),3) ≈ chebyshevutransform!(copy(X),Val(2),3) ≈ X̃ + + @test @inferred(chebyshevutransform(X)) ≈ @inferred(chebyshevutransform!(copy(X))) ≈ chebyshevutransform(chebyshevutransform(chebyshevutransform(X,1),2),3) + @test @inferred(chebyshevutransform(X,Val(2))) ≈ @inferred(chebyshevutransform!(copy(X),Val(2))) ≈ chebyshevutransform(chebyshevutransform(chebyshevutransform(X,Val(2),1),Val(2),2),Val(2),3) + end + + @testset "ichebyshevutransform" begin + for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = ichebyshevutransform(X[:,k,j]) end + @test @inferred(ichebyshevutransform(X,1)) ≈ @inferred(ichebyshevutransform!(copy(X),1)) ≈ X̃ + for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = ichebyshevutransform(X[k,:,j]) end + @test ichebyshevutransform(X,2) ≈ ichebyshevutransform!(copy(X),2) ≈ X̃ + for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = ichebyshevutransform(X[k,j,:]) end + @test ichebyshevutransform(X,3) ≈ ichebyshevutransform!(copy(X),3) ≈ X̃ + + for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = ichebyshevutransform(X[:,k,j],Val(2)) end + @test @inferred(ichebyshevutransform(X,Val(2),1)) ≈ @inferred(ichebyshevutransform!(copy(X),Val(2),1)) ≈ X̃ + for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = ichebyshevutransform(X[k,:,j],Val(2)) end + @test ichebyshevutransform(X,Val(2),2) ≈ ichebyshevutransform!(copy(X),Val(2),2) ≈ X̃ + for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = ichebyshevutransform(X[k,j,:],Val(2)) end + @test ichebyshevutransform(X,Val(2),3) ≈ ichebyshevutransform!(copy(X),Val(2),3) ≈ X̃ + + @test @inferred(ichebyshevutransform(X)) ≈ @inferred(ichebyshevutransform!(copy(X))) ≈ ichebyshevutransform(ichebyshevutransform(ichebyshevutransform(X,1),2),3) + @test @inferred(ichebyshevutransform(X,Val(2))) ≈ @inferred(ichebyshevutransform!(copy(X),Val(2))) ≈ ichebyshevutransform(ichebyshevutransform(ichebyshevutransform(X,Val(2),1),Val(2),2),Val(2),3) + + @test ichebyshevutransform(chebyshevutransform(X)) ≈ X + @test chebyshevutransform(ichebyshevutransform(X)) ≈ X + end + + X = randn(1,1,1) + @test chebyshevtransform!(copy(X), Val(1)) == ichebyshevtransform!(copy(X), Val(1)) == X + @test_throws ArgumentError chebyshevtransform!(copy(X), Val(2)) + @test_throws ArgumentError ichebyshevtransform!(copy(X), Val(2)) end - @testset "ichebyshevutransform" begin - for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = ichebyshevutransform(X[:,k,j]) end - @test @inferred(ichebyshevutransform(X,1)) ≈ @inferred(ichebyshevutransform!(copy(X),1)) ≈ X̃ - for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = ichebyshevutransform(X[k,:,j]) end - @test ichebyshevutransform(X,2) ≈ ichebyshevutransform!(copy(X),2) ≈ X̃ - for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = ichebyshevutransform(X[k,j,:]) end - @test ichebyshevutransform(X,3) ≈ ichebyshevutransform!(copy(X),3) ≈ X̃ - - for k = axes(X,2), j = axes(X,3) X̃[:,k,j] = ichebyshevutransform(X[:,k,j],Val(2)) end - @test @inferred(ichebyshevutransform(X,Val(2),1)) ≈ @inferred(ichebyshevutransform!(copy(X),Val(2),1)) ≈ X̃ - for k = axes(X,1), j = axes(X,3) X̃[k,:,j] = ichebyshevutransform(X[k,:,j],Val(2)) end - @test ichebyshevutransform(X,Val(2),2) ≈ ichebyshevutransform!(copy(X),Val(2),2) ≈ X̃ - for k = axes(X,1), j = axes(X,2) X̃[k,j,:] = ichebyshevutransform(X[k,j,:],Val(2)) end - @test ichebyshevutransform(X,Val(2),3) ≈ ichebyshevutransform!(copy(X),Val(2),3) ≈ X̃ - - @test @inferred(ichebyshevutransform(X)) ≈ @inferred(ichebyshevutransform!(copy(X))) ≈ ichebyshevutransform(ichebyshevutransform(ichebyshevutransform(X,1),2),3) - @test @inferred(ichebyshevutransform(X,Val(2))) ≈ @inferred(ichebyshevutransform!(copy(X),Val(2))) ≈ ichebyshevutransform(ichebyshevutransform(ichebyshevutransform(X,Val(2),1),Val(2),2),Val(2),3) - - @test ichebyshevutransform(chebyshevutransform(X)) ≈ X - @test chebyshevutransform(ichebyshevutransform(X)) ≈ X + @testset "4D" begin + X = randn(2,3,4,5) + X̃ = similar(X) + for trans in (chebyshevtransform, ichebyshevtransform, chebyshevutransform, ichebyshevutransform) + for k = axes(X,2), j = axes(X,3), l = axes(X,4) X̃[:,k,j,l] = trans(X[:,k,j,l]) end + @test @inferred(trans(X,1)) ≈ X̃ + @test @inferred(trans(X)) ≈ trans(trans(trans(trans(X,1),2),3),4) + end end - - X = randn(1,1,1) - @test chebyshevtransform!(copy(X), Val(1)) == ichebyshevtransform!(copy(X), Val(1)) == X - @test_throws ArgumentError chebyshevtransform!(copy(X), Val(2)) - @test_throws ArgumentError ichebyshevtransform!(copy(X), Val(2)) end - @testset "Integer" begin @test chebyshevtransform([1,2,3]) == chebyshevtransform([1.,2,3]) @test chebyshevtransform([1,2,3], Val(2)) == chebyshevtransform([1.,2,3], Val(2)) From a1024f208f64ee8a82b3c19ac4e5699c4a8e41d9 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 16 Nov 2023 10:47:32 +0000 Subject: [PATCH 12/25] Support planning Chebyshev via specifying type and size (#236) * Support planning Chebyshev via specifying type and size * Update Project.toml --- Project.toml | 2 +- src/chebyshevtransform.jl | 23 +++++++++++------------ test/chebyshevtests.jl | 6 ++++++ 3 files changed, 18 insertions(+), 13 deletions(-) diff --git a/Project.toml b/Project.toml index 3d45e5c3..9a985c5e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.15.13" +version = "0.15.14" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/src/chebyshevtransform.jl b/src/chebyshevtransform.jl index d55b52c3..5303868e 100644 --- a/src/chebyshevtransform.jl +++ b/src/chebyshevtransform.jl @@ -49,9 +49,6 @@ function plan_chebyshevtransform(x::AbstractArray{T,N}, ::Val{2}, dims...; kws.. ChebyshevTransformPlan{T,2}(FFTW.plan_r2r(x, SECONDKIND, dims...; kws...)) end -plan_chebyshevtransform!(x::AbstractArray, dims...; kws...) = plan_chebyshevtransform!(x, Val(1), dims...; kws...) -plan_chebyshevtransform(x::AbstractArray, dims...; kws...) = plan_chebyshevtransform(x, Val(1), dims...; kws...) - # convert x if necessary @inline _plan_mul!(y::AbstractArray{T}, P::Plan{T}, x::StridedArray{T}) where T = mul!(y, P, x) @@ -221,9 +218,6 @@ function plan_ichebyshevtransform(x::AbstractArray{T}, ::Val{2}, dims...; kws... inv(plan_chebyshevtransform(x, Val(2), dims...; kws...)) end -plan_ichebyshevtransform!(x::AbstractArray, dims...; kws...) = plan_ichebyshevtransform!(x, Val(1), dims...; kws...) -plan_ichebyshevtransform(x::AbstractArray, dims...; kws...) = plan_ichebyshevtransform(x, Val(1), dims...; kws...) - @inline function _icheb1_prescale!(d::Number, x::AbstractArray) lmul_dim_begin!(2, d, x) x @@ -369,9 +363,6 @@ function plan_chebyshevutransform(x::AbstractArray{T,N}, ::Val{2}, dims...; kws. ChebyshevUTransformPlan{T,2}(FFTW.plan_r2r(x, USECONDKIND, dims...; kws...)) end -plan_chebyshevutransform!(x::AbstractArray, dims...; kws...) = plan_chebyshevutransform!(x, Val(1), dims...; kws...) -plan_chebyshevutransform(x::AbstractArray, dims...; kws...) = plan_chebyshevutransform(x, Val(1), dims...; kws...) - _permfirst(d, N) = [d; 1:d-1; d+1:N] @@ -526,9 +517,6 @@ function plan_ichebyshevutransform(x::AbstractArray{T,N}, ::Val{2}, dims...; kws end -plan_ichebyshevutransform!(x::AbstractArray, dims...; kws...) = plan_ichebyshevutransform!(x, Val(1), dims...; kws...) -plan_ichebyshevutransform(x::AbstractArray, dims...; kws...) = plan_ichebyshevutransform(x, Val(1), dims...; kws...) - # second kind Chebyshev transforms share a plan with their inverse # so we support this via inv inv(P::ChebyshevUTransformPlan{T,2}) where {T} = IChebyshevUTransformPlan{T,2}(P.plan) @@ -744,3 +732,14 @@ end copyto!(x, IChebyshevTransformPlan{T,1,Nothing,false,N,R}() * x) # *(P::IChebyshevTransformPlan{T,SECONDKIND,false,Nothing}, x::AbstractVector{T}) where T = # IChebyshevTransformPlan{T,SECONDKIND,true,Nothing}() * copy(x) + + +for pln in (:plan_chebyshevtransform!, :plan_chebyshevtransform, + :plan_chebyshevutransform!, :plan_chebyshevutransform, + :plan_ichebyshevutransform, :plan_ichebyshevutransform!, + :plan_ichebyshevtransform, :plan_ichebyshevtransform!) + @eval begin + $pln(x::AbstractArray, dims...; kws...) = $pln(x, Val(1), dims...; kws...) + $pln(::Type{T}, szs, dims...; kwds...) where T = $pln(Array{T}(undef, szs...), dims...; kwds...) + end +end \ No newline at end of file diff --git a/test/chebyshevtests.jl b/test/chebyshevtests.jl index 2a82978c..3112472d 100644 --- a/test/chebyshevtests.jl +++ b/test/chebyshevtests.jl @@ -478,4 +478,10 @@ using FastTransforms, Test @test_throws ErrorException plan_chebyshevtransform(randn(5)) * randn(5,5) @test_throws ErrorException plan_ichebyshevtransform(randn(5)) * randn(5,5) end + + @testset "plan via size" begin + X = randn(3,4) + p = plan_chebyshevtransform(Float64, (3,4)) + @test p * X == chebyshevtransform(X) + end end From e3d47a2d429b2c4425a520fb4bb4e3da4e3f10d4 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Thu, 16 Nov 2023 14:40:41 +0000 Subject: [PATCH 13/25] Fix bug in ultra2ultra with integer difference (#237) --- Project.toml | 2 +- src/toeplitzhankel.jl | 2 +- test/toeplitzhankeltests.jl | 1 + 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 9a985c5e..022e0477 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.15.14" +version = "0.15.15" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/src/toeplitzhankel.jl b/src/toeplitzhankel.jl index 383ff9ae..f7825466 100644 --- a/src/toeplitzhankel.jl +++ b/src/toeplitzhankel.jl @@ -299,7 +299,7 @@ function plan_th_ultra2ultra!(::Type{S}, mn, λ₁, λ₂, dims) where {S} if isapproxinteger(λ₂ - λ₁) # TODO: don't make extra plan - plans = typeof(_good_plan_th_ultra2ultra!(S, mn, λ₁+0.1, λ₂, dims))[] + plans = typeof(_good_plan_th_ultra2ultra!(S, mn, λ₂+0.1, λ₂, dims))[] else plans = [_good_plan_th_ultra2ultra!(S, mn, λ₁, c, dims)] end diff --git a/test/toeplitzhankeltests.jl b/test/toeplitzhankeltests.jl index b90ce94b..a8ebdc40 100644 --- a/test/toeplitzhankeltests.jl +++ b/test/toeplitzhankeltests.jl @@ -15,6 +15,7 @@ Random.seed!(0) @test th_ultra2ultra(x,1, 2) ≈ lib_ultra2ultra(x, 1, 2) @test th_ultra2ultra(x,0.1, 2.2) ≈ lib_ultra2ultra(x, 0.1, 2.2) @test th_ultra2ultra(x, 2.2, 0.1) ≈ lib_ultra2ultra(x, 2.2, 0.1) + @test th_ultra2ultra(x, 1, 3) ≈ lib_ultra2ultra(x, 1, 3) @test @inferred(th_jac2jac(x,0.1, 0.2,0.1,0.4)) ≈ lib_jac2jac(x, 0.1, 0.2,0.1,0.4) @test th_jac2jac(x,0.1, 0.2,0.3,0.2) ≈ lib_jac2jac(x, 0.1, 0.2,0.3,0.2) @test th_jac2jac(x,0.1, 0.2,0.3,0.4) ≈ lib_jac2jac(x, 0.1, 0.2,0.3,0.4) From 9f135a895054c81fc105fd10058bbc8cfca1467d Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Mon, 4 Dec 2023 16:55:19 +0000 Subject: [PATCH 14/25] Support size and \ for cheb2leg (#239) --- Project.toml | 2 +- src/toeplitzhankel.jl | 8 ++++++++ test/toeplitzhankeltests.jl | 19 +++++++++++++++++++ 3 files changed, 28 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 022e0477..2f5b3613 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.15.15" +version = "0.15.16" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/src/toeplitzhankel.jl b/src/toeplitzhankel.jl index f7825466..78cc43fa 100644 --- a/src/toeplitzhankel.jl +++ b/src/toeplitzhankel.jl @@ -33,6 +33,8 @@ ToeplitzHankelPlan{S,N,M}(T::TP, L::LowR, R::LowR, dims::Dims) where {S,N,M,LowR ToeplitzHankelPlan{S,N}(T, L, R, dims) where {S,N} = ToeplitzHankelPlan{S,N,N+1}(T, L, R, dims) ToeplitzHankelPlan(T::ToeplitzPlan{S,M}, L::Matrix, R::Matrix, dims=1) where {S,M} = ToeplitzHankelPlan{S,M-1,M}((T,), (L,), (R,), dims) +size(TH::ToeplitzHankelPlan) = size(first(TH.T)) + _reshape_broadcast(d, R, ::Val{N}, M) where N = reshape(R,ntuple(k -> k == d ? size(R,1) : 1, Val(N))...,M) function _th_applymul!(d, v::AbstractArray{<:Any,N}, T, L, R, tmp) where N @@ -52,6 +54,8 @@ function *(P::ToeplitzHankelPlan{<:Any,N}, v::AbstractArray{<:Any,N}) where N v end +*(P::ToeplitzHankelPlan, v::AbstractArray) = error("plan applied to wrong-sized array") + # partial cholesky for a Hankel matrix @@ -158,6 +162,10 @@ function *(P::ChebyshevToLegendrePlanTH, V::AbstractArray{<:Any,N}) where N V end +_add1tod(d::Integer, a, b...) = d == 1 ? (a+1, b...) : (a, _add1tod(d-1, b...)...) +_add1tod(d, a, b...) = _add1tod(first(d), a, b...) +size(P::ChebyshevToLegendrePlanTH) = Base.front(_add1tod(P.toeplitzhankel.dims, size(first(P.toeplitzhankel.T))...)) +inv(P::ChebyshevToLegendrePlanTH{T}) where T = plan_th_leg2cheb!(T, size(P), P.toeplitzhankel.dims) function _leg2chebTH_TLC(::Type{S}, mn, d) where S diff --git a/test/toeplitzhankeltests.jl b/test/toeplitzhankeltests.jl index a8ebdc40..0b8731bf 100644 --- a/test/toeplitzhankeltests.jl +++ b/test/toeplitzhankeltests.jl @@ -164,4 +164,23 @@ Random.seed!(0) @test M ≈ Y end end + + @testset "inv" begin + x = randn(10) + pl = plan_th_cheb2leg!(x) + @test size(pl) == (10,) + @test pl\(pl*x) ≈ x + + X = randn(10,3) + for pl in (plan_th_cheb2leg!(X), plan_th_cheb2leg!(X, 1), plan_th_cheb2leg!(X, 2)) + @test size(pl) == (10,3) + @test pl\(pl*copy(X)) ≈ X + end + + X = randn(10,3,5) + for pl in (plan_th_cheb2leg!(X), plan_th_cheb2leg!(X, 1), plan_th_cheb2leg!(X, 2), plan_th_cheb2leg!(X, 3)) + @test size(pl) == (10,3,5) + @test pl\(pl*copy(X)) ≈ X + end + end end \ No newline at end of file From 7ca43ceb57c1c061f9c2fea635a04cd7ced7254b Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Fri, 15 Mar 2024 10:43:23 -0500 Subject: [PATCH 15/25] add Symmetric(Banded)ToeplitzPlusHankel --- Project.toml | 4 +- src/FastTransforms.jl | 8 +- src/SymmetricToeplitzPlusHankel.jl | 135 +++++++++++++++++++++++ test/runtests.jl | 3 +- test/symmetrictoeplitzplushankeltests.jl | 39 +++++++ 5 files changed, 185 insertions(+), 4 deletions(-) create mode 100644 src/SymmetricToeplitzPlusHankel.jl create mode 100644 test/symmetrictoeplitzplushankeltests.jl diff --git a/Project.toml b/Project.toml index 2f5b3613..8bba0be1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,9 +1,10 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.15.16" +version = "0.16.0" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" +BandedMatrices = "aae01518-5342-5314-be14-df237901396f" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" FastTransforms_jll = "34b6f7d7-08f9-5794-9e10-3819e4c7e49a" @@ -17,6 +18,7 @@ ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24" [compat] AbstractFFTs = "1.0" +BandedMatrices = "1.5" FFTW = "1.7" FastGaussQuadrature = "0.4, 0.5, 1" FastTransforms_jll = "0.6.2" diff --git a/src/FastTransforms.jl b/src/FastTransforms.jl index e9bf6fac..82c34990 100644 --- a/src/FastTransforms.jl +++ b/src/FastTransforms.jl @@ -1,6 +1,6 @@ module FastTransforms -using FastGaussQuadrature, FillArrays, LinearAlgebra, +using BandedMatrices, FastGaussQuadrature, FillArrays, LinearAlgebra, Reexport, SpecialFunctions, ToeplitzMatrices @reexport using AbstractFFTs @@ -19,6 +19,8 @@ import AbstractFFTs: Plan, ScaledPlan, fftshift, ifftshift, rfft_output_size, brfft_output_size, normalization +import BandedMatrices: bandwidths + import FFTW: dct, dct!, idct, idct!, plan_dct!, plan_idct!, plan_dct, plan_idct, fftwNumber @@ -26,7 +28,7 @@ import FastGaussQuadrature: unweightedgausshermite import FillArrays: AbstractFill, getindex_value -import LinearAlgebra: mul!, lmul!, ldiv! +import LinearAlgebra: mul!, lmul!, ldiv!, cholesky import GenericFFT: interlace # imported in downstream packages @@ -112,6 +114,8 @@ include("specialfunctions.jl") include("toeplitzplans.jl") include("toeplitzhankel.jl") +include("SymmetricToeplitzPlusHankel.jl") + # following use libfasttransforms by default for f in (:jac2jac, :lag2lag, :jac2ultra, :ultra2jac, :jac2cheb, diff --git a/src/SymmetricToeplitzPlusHankel.jl b/src/SymmetricToeplitzPlusHankel.jl new file mode 100644 index 00000000..a9bc4013 --- /dev/null +++ b/src/SymmetricToeplitzPlusHankel.jl @@ -0,0 +1,135 @@ +struct SymmetricToeplitzPlusHankel{T} <: AbstractMatrix{T} + v::Vector{T} + n::Int +end + +function SymmetricToeplitzPlusHankel(v::Vector{T}) where T + n = (length(v)+1)÷2 + SymmetricToeplitzPlusHankel{T}(v, n) +end + +size(A::SymmetricToeplitzPlusHankel{T}) where T = (A.n, A.n) +getindex(A::SymmetricToeplitzPlusHankel{T}, i::Integer, j::Integer) where T = A.v[abs(i-j)+1] + A.v[i+j-1] + +struct SymmetricBandedToeplitzPlusHankel{T} <: BandedMatrices.AbstractBandedMatrix{T} + v::Vector{T} + n::Int + b::Int +end + +function SymmetricBandedToeplitzPlusHankel(v::Vector{T}, n::Integer) where T + SymmetricBandedToeplitzPlusHankel{T}(v, n, length(v)-1) +end + +size(A::SymmetricBandedToeplitzPlusHankel{T}) where T = (A.n, A.n) +function getindex(A::SymmetricBandedToeplitzPlusHankel{T}, i::Integer, j::Integer) where T + v = A.v + if abs(i-j) < length(v) + if i+j-1 ≤ length(v) + v[abs(i-j)+1] + v[i+j-1] + else + v[abs(i-j)+1] + end + else + zero(T) + end +end +bandwidths(A::SymmetricBandedToeplitzPlusHankel{T}) where T = (A.b, A.b) + +# +# Jac*W-W*Jac' = G*J*G' +# This returns G and J, where J = [0 I; -I 0], respecting the skew-symmetry of the right-hand side. +# +function compute_skew_generators(A::SymmetricToeplitzPlusHankel{T}) where T + v = A.v + n = size(A, 1) + J = [zero(T) one(T); -one(T) zero(T)] + G = zeros(T, n, 2) + G[n, 1] = one(T) + u2 = reverse(v[2:n+1]) + u2[1:n-1] .+= v[n+1:2n-1] + G[:, 2] .= -u2 + G, J +end + +function cholesky(A::SymmetricToeplitzPlusHankel{T}) where T + n = size(A, 1) + G, J = compute_skew_generators(A) + L = zeros(T, n, n) + r = A[:, 1] + r2 = zeros(T, n) + l = zeros(T, n) + v = zeros(T, n) + col1 = zeros(T, n) + STPHcholesky!(L, G, r, r2, l, v, col1, n) + return Cholesky(L, 'L', 0) +end + +function STPHcholesky!(L::Matrix{T}, G, r, r2, l, v, col1, n) where T + @inbounds @simd for k in 1:n-1 + x = sqrt(r[1]) + for j in 1:n-k+1 + L[j+k-1, k] = l[j] = r[j]/x + end + for j in 1:n-k+1 + v[j] = G[j, 1]*G[1,2]-G[j,2]*G[1,1] + end + F12 = k == 1 ? T(2) : T(1) + r2[1] = (r[2] - v[1])/F12 + for j in 2:n-k + r2[j] = (r[j+1]+r[j-1] + r[1]*col1[j] - col1[1]*r[j] - v[j])/F12 + end + r2[n-k+1] = (r[n-k] + r[1]*col1[n-k+1] - col1[1]*r[n-k+1] - v[n-k+1])/F12 + cst = r[2]/x + for j in 1:n-k + r[j] = r2[j+1] - cst*l[j+1] + end + for j in 1:n-k + col1[j] = -F12/x*l[j+1] + end + c1 = G[1, 1] + c2 = G[1, 2] + for j in 1:n-k + G[j, 1] = G[j+1, 1] - l[j+1]*c1/x + G[j, 2] = G[j+1, 2] - l[j+1]*c2/x + end + end + L[n, n] = sqrt(r[1]) +end + +function cholesky(A::SymmetricBandedToeplitzPlusHankel{T}) where T + n = size(A, 1) + b = A.b + R = BandedMatrix{T}(undef, (n, n), (0, bandwidth(A, 2))) + r = A[1:b+2, 1] + r2 = zeros(T, b+3) + l = zeros(T, b+3) + col1 = zeros(T, b+2) + SBTPHcholesky!(R, r, r2, l, col1, n, b) + return Cholesky(R, 'U', 0) +end + +function SBTPHcholesky!(R::BandedMatrix{T}, r, r2, l, col1, n, b) where T + @inbounds @simd for k in 1:n + x = sqrt(r[1]) + for j in 1:b+1 + l[j] = r[j]/x + end + for j in 1:min(n-k+1, b+1) + R[k, j+k-1] = l[j] + end + F12 = k == 1 ? T(2) : T(1) + r2[1] = r[2]/F12 + for j in 2:b+1 + r2[j] = (r[j+1]+r[j-1] + r[1]*col1[j] - col1[1]*r[j])/F12 + end + r2[b+2] = (r[b+1] + r[1]*col1[b+2] - col1[1]*r[b+2])/F12 + cst = r[2]/x + for j in 1:b+2 + r[j] = r2[j+1] - cst*l[j+1] + end + for j in 1:b+2 + col1[j] = -F12/x*l[j+1] + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 1efb21e9..de16f36b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,4 +10,5 @@ include("gaunttests.jl") include("hermitetests.jl") include("clenshawtests.jl") include("toeplitzplanstests.jl") -include("toeplitzhankeltests.jl") \ No newline at end of file +include("toeplitzhankeltests.jl") +include("symmetrictoeplitzplushankeltests.jl") diff --git a/test/symmetrictoeplitzplushankeltests.jl b/test/symmetrictoeplitzplushankeltests.jl new file mode 100644 index 00000000..c26c38d3 --- /dev/null +++ b/test/symmetrictoeplitzplushankeltests.jl @@ -0,0 +1,39 @@ +using BandedMatrices, FastTransforms, LinearAlgebra, ToeplitzMatrices, Test + +import FastTransforms: SymmetricToeplitzPlusHankel, SymmetricBandedToeplitzPlusHankel + +@testset "SymmetricToeplitzPlusHankel" begin + n = 128 + for T in (Float32, Float64, BigFloat) + μ = -FastTransforms.chebyshevlogmoments1(T, 2n-1) + μ[1] += 1 + W = SymmetricToeplitzPlusHankel(μ/2) + SMW = Symmetric(Matrix(W)) + @test W ≈ SymmetricToeplitz(μ[1:(length(μ)+1)÷2]/2) + Hankel(μ/2) + L = cholesky(W).L + R = cholesky(SMW).U + @test L*L' ≈ W + @test L' ≈ R + end +end + +@testset "SymmetricBandedToeplitzPlusHankel" begin + n = 1024 + for T in (Float32, Float64) + μ = T[1.875; 0.00390625; 0.5; 0.0009765625; 0.0625] + W = SymmetricBandedToeplitzPlusHankel(μ/2, n) + SBW = Symmetric(BandedMatrix(W)) + W1 = SymmetricToeplitzPlusHankel([μ/2; zeros(2n-1-length(μ))]) + SMW = Symmetric(Matrix(W)) + U = cholesky(SMW).U + L = cholesky(W1).L + UB = cholesky(SBW).U + R = cholesky(W).U + @test L*L' ≈ W + @test UB'UB ≈ W + @test R'R ≈ W + @test UB ≈ U + @test L' ≈ U + @test R ≈ U + end +end From 324c2230c673ad270291fb2731e714ce9317757a Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Fri, 15 Mar 2024 15:05:15 -0500 Subject: [PATCH 16/25] add absolute value moments --- Project.toml | 2 +- src/specialfunctions.jl | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 8bba0be1..97d36acb 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "FastTransforms" uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c" -version = "0.16.0" +version = "0.16.1" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/src/specialfunctions.jl b/src/specialfunctions.jl index 66ef53c2..b5915ab7 100644 --- a/src/specialfunctions.jl +++ b/src/specialfunctions.jl @@ -307,6 +307,23 @@ function chebyshevlogmoments1(::Type{T}, N::Int) where T μ end +""" +Modified Chebyshev moments of the first kind with respect to the absolute value weight: + +```math + \\int_{-1}^{+1} T_n(x) |x|{\\rm\\,d}x. +``` +""" +function chebyshevabsmoments1(::Type{T}, N::Int) where T + μ = zeros(T, N) + if N > 0 + for i=0:4:N-1 + @inbounds μ[i+1] = -T(1)/T((i÷2)^2-1) + end + end + μ +end + """ Modified Chebyshev moments of the second kind: @@ -359,6 +376,23 @@ function chebyshevlogmoments2(::Type{T}, N::Int) where T μ end +""" +Modified Chebyshev moments of the second kind with respect to the absolute value weight: + +```math + \\int_{-1}^{+1} U_n(x) |x|{\\rm\\,d}x. +``` +""" +function chebyshevabsmoments2(::Type{T}, N::Int) where T + μ = chebyshevabsmoments1(T, N) + if N > 1 + μ[2] *= two(T) + for i=1:N-2 + @inbounds μ[i+2] = 2μ[i+2] + μ[i] + end + end + μ +end function sphrand(::Type{T}, m::Int, n::Int) where T A = zeros(T, m, n) From 2c763590b22f8c38e3a60c0e296776c66f49e411 Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Thu, 21 Mar 2024 11:04:36 -0500 Subject: [PATCH 17/25] create halfrange example --- docs/Project.toml | 1 + docs/make.jl | 2 ++ examples/halfrange.jl | 62 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 65 insertions(+) create mode 100644 examples/halfrange.jl diff --git a/docs/Project.toml b/docs/Project.toml index cabe3c60..c96b3cd3 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,7 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c" +LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" diff --git a/docs/make.jl b/docs/make.jl index 77842971..dd530e8f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -10,6 +10,7 @@ examples = [ "automaticdifferentiation.jl", "chebyshev.jl", "disk.jl", + "halfrange.jl", "nonlocaldiffusion.jl", "padua.jl", "sphere.jl", @@ -43,6 +44,7 @@ makedocs( "generated/automaticdifferentiation.md", "generated/chebyshev.md", "generated/disk.md", + "generated/halfrange.md", "generated/nonlocaldiffusion.md", "generated/padua.md", "generated/sphere.md", diff --git a/examples/halfrange.jl b/examples/halfrange.jl new file mode 100644 index 00000000..6d9f6a69 --- /dev/null +++ b/examples/halfrange.jl @@ -0,0 +1,62 @@ +# # Half-range Chebyshev polynomials +# In [this paper](https://doi.org/10.1137/090752456), [Daan Huybrechs](https://github.com/daanhb) introduced the so-called half-range Chebyshev polynomials +# as the non-classical orthogonal polynomials with respect to the inner product: +# ```math +# \langle f, g \rangle = \int_0^1 f(x) g(x)\frac{{\rm d} x}{\sqrt{1-x^2}}. +# ``` +# By the variable transformation $y = 2x-1$, the resulting polynomials can be related to +# orthogonal polynomials on $(-1,1)$ with the Jacobi weight $(1-y)^{-\frac{1}{2}}$ modified by the weight $(3+y)^{-\frac{1}{2}}$. +# +# We shall use the fact that: +# ```math +# \frac{1}{\sqrt{3+y}} = \sqrt{\frac{2}{3+\sqrt{8}}}\sum_{n=0}^\infty P_n(y) \left(\frac{-1}{3+\sqrt{8}}\right)^n, +# ``` +# and results from [this paper](https://arxiv.org/abs/2302.08448) to consider the half-range Chebyshev polynomials as +# modifications of the Jacobi polynomials $P_n^{(-\frac{1}{2},0)}(y)$. + +using FastTransforms, LinearAlgebra, Plots, LaTeXStrings +const GENFIGS = joinpath(pkgdir(FastTransforms), "docs/src/generated") +!isdir(GENFIGS) && mkdir(GENFIGS) +plotlyjs() + +# We truncate the generating function to ensure an absolute error of `eps()`: +z = -1/(3+sqrt(8)) +K = sqrt(-2z) +N = log(abs(z), eps()*(1-abs(z))/K) - 1 +d = K .* z .^(0:N) + +# Then, we convert this representation to the expansion in Jacobi polynomials $P_n^{(-\frac{1}{2}, 0)}(y)$: +u = jac2jac(d, 0.0, 0.0, -0.5, 0.0; norm1 = false, norm2 = true) + +# Our working polynomial degree will be: +n = 100 + +# We compute the connection coefficients between the modified orthogonal polynomials and the Jacobi polynomials: +P = plan_modifiedjac2jac(Float64, n+1, -0.5, 0.0, u) + +# We store the connection to first kind Chebyshev polynomials: +P1 = plan_jac2cheb(Float64, n+1, -0.5, 0.0; normjac = true) + +# We compute the Chebyshev series for the degree-$k\le n$ modified polynomial and its values at the Chebyshev points: +q = k -> lmul!(P1, lmul!(P, [zeros(k); 1.0; zeros(n-k)])) +qvals = k-> ichebyshevtransform(q(k)) + +# With the symmetric Jacobi matrix for $P_n^{(-\frac{1}{2}, 0)}(y)$ and the modified plan, we may compute the modified Jacobi matrix and the corresponding roots (as eigenvalues): +XP = SymTridiagonal([-inv((4n-1)*(4n-5)) for n in 1:n+1], [4n*(2n-1)/(4n-1)/sqrt((4n-3)*(4n+1)) for n in 1:n]) +XQ = FastTransforms.modified_jacobi_matrix(P, XP) + +# And we plot: +x = (chebyshevpoints(Float64, n+1, Val(1)) .+ 1 ) ./ 2 +p = plot(x, qvals(0); linewidth=2.0, legend = false, xlim=(0,1), xlabel=L"x", + ylabel=L"T^h_n(x)", title="Half-Range Chebyshev Polynomials and Their Roots", + extra_plot_kwargs = KW(:include_mathjax => "cdn")) +for k in 1:10 + λ = (eigvals(SymTridiagonal(XQ.dv[1:k], XQ.ev[1:k-1])) .+ 1) ./ 2 + plot!(x, qvals(k); linewidth=2.0, color=palette(:default)[k+1]) + scatter!(λ, zero(λ); markersize=2.5, color=palette(:default)[k+1]) +end +p +#savefig(joinpath(GENFIGS, "halfrange.html")) +###```@raw html +### +###``` From 1bfa33b8856c8a77822dc74fe448deaeaa39d68e Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Thu, 21 Mar 2024 11:05:55 -0500 Subject: [PATCH 18/25] uncomment the savefig --- examples/halfrange.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/halfrange.jl b/examples/halfrange.jl index 6d9f6a69..eec945ab 100644 --- a/examples/halfrange.jl +++ b/examples/halfrange.jl @@ -56,7 +56,7 @@ for k in 1:10 scatter!(λ, zero(λ); markersize=2.5, color=palette(:default)[k+1]) end p -#savefig(joinpath(GENFIGS, "halfrange.html")) +savefig(joinpath(GENFIGS, "halfrange.html")) ###```@raw html ### ###``` From 04d8e8718ee41024fd47dd2aeddd028bc49f586b Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Thu, 21 Mar 2024 11:19:28 -0500 Subject: [PATCH 19/25] curtail large output of SymTridiagonal --- examples/halfrange.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/halfrange.jl b/examples/halfrange.jl index eec945ab..7205b237 100644 --- a/examples/halfrange.jl +++ b/examples/halfrange.jl @@ -43,7 +43,7 @@ qvals = k-> ichebyshevtransform(q(k)) # With the symmetric Jacobi matrix for $P_n^{(-\frac{1}{2}, 0)}(y)$ and the modified plan, we may compute the modified Jacobi matrix and the corresponding roots (as eigenvalues): XP = SymTridiagonal([-inv((4n-1)*(4n-5)) for n in 1:n+1], [4n*(2n-1)/(4n-1)/sqrt((4n-3)*(4n+1)) for n in 1:n]) -XQ = FastTransforms.modified_jacobi_matrix(P, XP) +XQ = FastTransforms.modified_jacobi_matrix(P, XP); # And we plot: x = (chebyshevpoints(Float64, n+1, Val(1)) .+ 1 ) ./ 2 From 78c98fb664f60d9d9053a98db293cfd99567677b Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Thu, 21 Mar 2024 12:39:33 -0500 Subject: [PATCH 20/25] include banded diff --- examples/halfrange.jl | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/examples/halfrange.jl b/examples/halfrange.jl index 7205b237..5fce0ea3 100644 --- a/examples/halfrange.jl +++ b/examples/halfrange.jl @@ -1,6 +1,6 @@ # # Half-range Chebyshev polynomials # In [this paper](https://doi.org/10.1137/090752456), [Daan Huybrechs](https://github.com/daanhb) introduced the so-called half-range Chebyshev polynomials -# as the non-classical orthogonal polynomials with respect to the inner product: +# as the semi-classical orthogonal polynomials with respect to the inner product: # ```math # \langle f, g \rangle = \int_0^1 f(x) g(x)\frac{{\rm d} x}{\sqrt{1-x^2}}. # ``` @@ -19,10 +19,10 @@ const GENFIGS = joinpath(pkgdir(FastTransforms), "docs/src/generated") !isdir(GENFIGS) && mkdir(GENFIGS) plotlyjs() -# We truncate the generating function to ensure an absolute error of `eps()`: +# We truncate the generating function to ensure a relative error less than `eps()` in the uniform norm on $(-1,1)$: z = -1/(3+sqrt(8)) K = sqrt(-2z) -N = log(abs(z), eps()*(1-abs(z))/K) - 1 +N = ceil(Int, log(abs(z), eps()/2*(1-abs(z))/K) - 1) d = K .* z .^(0:N) # Then, we convert this representation to the expansion in Jacobi polynomials $P_n^{(-\frac{1}{2}, 0)}(y)$: @@ -60,3 +60,16 @@ savefig(joinpath(GENFIGS, "halfrange.html")) ###```@raw html ### ###``` + +# By [Theorem 2.20](https://arxiv.org/abs/2302.08448) it turns out that the *derivatives* of the half-range Chebyshev polynomials are a linear combination of at most two polynomials orthogonal with respect to $\sqrt{(3+y)(1-y)}(1+y)$ on $(-1,1)$. This fact enables us to compute the banded differentiation matrix: +v̂ = 3*u+XP[1:N+1,1:N]*u +v = jac2jac(v̂, -0.5, 0.0, 0.5, 1.0; norm1 = true, norm2 = true) +function threshold!(A::AbstractArray, ϵ) + for i in eachindex(A) + if abs(A[i]) < ϵ A[i] = 0 end + end + A +end +P′ = plan_modifiedjac2jac(Float64, n+1, 0.5, 1.0, v) +DP = UpperTriangular(diagm(1=>[sqrt(n*(n+1/2)) for n in 1:n])) # The classical differentiation matrix representing 𝒟 P^{(-1/2,0)}(y) = P^{(1/2,1)}(y) D_P. +DQ = UpperTriangular(threshold!(P′\(DP*(P*I)), 100eps())) # The semi-classical differentiation matrix representing 𝒟 Q(y) = Q̂(y) D_Q. From f0b71c6a9cb746234870905d07a50263b055fe8d Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Thu, 21 Mar 2024 12:50:57 -0500 Subject: [PATCH 21/25] add coefficients of the right sizes fix output sizes --- examples/halfrange.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/halfrange.jl b/examples/halfrange.jl index 5fce0ea3..f5152742 100644 --- a/examples/halfrange.jl +++ b/examples/halfrange.jl @@ -43,7 +43,8 @@ qvals = k-> ichebyshevtransform(q(k)) # With the symmetric Jacobi matrix for $P_n^{(-\frac{1}{2}, 0)}(y)$ and the modified plan, we may compute the modified Jacobi matrix and the corresponding roots (as eigenvalues): XP = SymTridiagonal([-inv((4n-1)*(4n-5)) for n in 1:n+1], [4n*(2n-1)/(4n-1)/sqrt((4n-3)*(4n+1)) for n in 1:n]) -XQ = FastTransforms.modified_jacobi_matrix(P, XP); +XQ = FastTransforms.modified_jacobi_matrix(P, XP) +SymTridiagonal(XQ.dv[1:10], XQ.ev[1:9]) # And we plot: x = (chebyshevpoints(Float64, n+1, Val(1)) .+ 1 ) ./ 2 @@ -62,7 +63,7 @@ savefig(joinpath(GENFIGS, "halfrange.html")) ###``` # By [Theorem 2.20](https://arxiv.org/abs/2302.08448) it turns out that the *derivatives* of the half-range Chebyshev polynomials are a linear combination of at most two polynomials orthogonal with respect to $\sqrt{(3+y)(1-y)}(1+y)$ on $(-1,1)$. This fact enables us to compute the banded differentiation matrix: -v̂ = 3*u+XP[1:N+1,1:N]*u +v̂ = 3*[u; 0]+XP[1:N+1, 1:N]*u v = jac2jac(v̂, -0.5, 0.0, 0.5, 1.0; norm1 = true, norm2 = true) function threshold!(A::AbstractArray, ϵ) for i in eachindex(A) @@ -73,3 +74,4 @@ end P′ = plan_modifiedjac2jac(Float64, n+1, 0.5, 1.0, v) DP = UpperTriangular(diagm(1=>[sqrt(n*(n+1/2)) for n in 1:n])) # The classical differentiation matrix representing 𝒟 P^{(-1/2,0)}(y) = P^{(1/2,1)}(y) D_P. DQ = UpperTriangular(threshold!(P′\(DP*(P*I)), 100eps())) # The semi-classical differentiation matrix representing 𝒟 Q(y) = Q̂(y) D_Q. +UpperTriangular(DQ[1:10,1:10]) From b2cd0e813c6730bbfec66c3b37b4340daf261fca Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Thu, 21 Mar 2024 13:04:58 -0500 Subject: [PATCH 22/25] update indexing --- examples/halfrange.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/halfrange.jl b/examples/halfrange.jl index f5152742..e9de076a 100644 --- a/examples/halfrange.jl +++ b/examples/halfrange.jl @@ -63,7 +63,7 @@ savefig(joinpath(GENFIGS, "halfrange.html")) ###``` # By [Theorem 2.20](https://arxiv.org/abs/2302.08448) it turns out that the *derivatives* of the half-range Chebyshev polynomials are a linear combination of at most two polynomials orthogonal with respect to $\sqrt{(3+y)(1-y)}(1+y)$ on $(-1,1)$. This fact enables us to compute the banded differentiation matrix: -v̂ = 3*[u; 0]+XP[1:N+1, 1:N]*u +v̂ = 3*[u; 0]+XP[1:N+2, 1:N+1]*u v = jac2jac(v̂, -0.5, 0.0, 0.5, 1.0; norm1 = true, norm2 = true) function threshold!(A::AbstractArray, ϵ) for i in eachindex(A) From 9c379b92b6653a9e50ae38ee887073a749f0b555 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sun, 28 Apr 2024 16:03:16 +0530 Subject: [PATCH 23/25] Fix type-promotion in BigInt transforms (#244) --- src/chebyshevtransform.jl | 4 ++-- test/chebyshevtests.jl | 4 ++++ 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/src/chebyshevtransform.jl b/src/chebyshevtransform.jl index 5303868e..7920e6d0 100644 --- a/src/chebyshevtransform.jl +++ b/src/chebyshevtransform.jl @@ -598,7 +598,7 @@ ichebyshevutransform(x, dims...; kwds...) = plan_ichebyshevutransform(x, dims... ## Code generation for integer inputs for func in (:chebyshevtransform,:ichebyshevtransform,:chebyshevutransform,:ichebyshevutransform) - @eval $func(x::AbstractVector{T}, dims...; kwds...) where {T<:Integer} = $func(convert(AbstractVector{Float64},x), dims...; kwds...) + @eval $func(x::AbstractVector{T}, dims...; kwds...) where {T<:Integer} = $func(convert(AbstractVector{float(T)},x), dims...; kwds...) end @@ -742,4 +742,4 @@ for pln in (:plan_chebyshevtransform!, :plan_chebyshevtransform, $pln(x::AbstractArray, dims...; kws...) = $pln(x, Val(1), dims...; kws...) $pln(::Type{T}, szs, dims...; kwds...) where T = $pln(Array{T}(undef, szs...), dims...; kwds...) end -end \ No newline at end of file +end diff --git a/test/chebyshevtests.jl b/test/chebyshevtests.jl index 3112472d..614f9c6d 100644 --- a/test/chebyshevtests.jl +++ b/test/chebyshevtests.jl @@ -443,6 +443,10 @@ using FastTransforms, Test @test plan_chebyshevtransform!(x)copy(x) ≈ chebyshevtransform(x) @test plan_ichebyshevtransform!(x)copy(x) ≈ ichebyshevtransform(x) end + @testset "BigInt" begin + x = big(10)^400 .+ BigInt[1,2,3] + @test ichebyshevtransform(chebyshevtransform(x)) ≈ x + end @testset "immutable vectors" begin F = plan_chebyshevtransform([1.,2,3]) From 8caad8184c15c2e26d27643b864d56d96a92a899 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 7 May 2024 12:35:32 +0530 Subject: [PATCH 24/25] Add PkgEval badge to readme (#245) --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index ac08bd9e..0686794d 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,7 @@ # FastTransforms.jl [![Build Status](https://github.com/JuliaApproximation/FastTransforms.jl/workflows/CI/badge.svg)](https://github.com/JuliaApproximation/FastTransforms.jl/actions?query=workflow%3ACI) [![codecov](https://codecov.io/gh/JuliaApproximation/FastTransforms.jl/branch/master/graph/badge.svg?token=BxTvSNgmLL)](https://codecov.io/gh/JuliaApproximation/FastTransforms.jl) [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaApproximation.github.io/FastTransforms.jl/stable) [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaApproximation.github.io/FastTransforms.jl/dev) +[![pkgeval](https://juliahub.com/docs/General/FastTransforms/stable/pkgeval.svg)](https://juliaci.github.io/NanosoldierReports/pkgeval_badges/report.html) `FastTransforms.jl` allows the user to conveniently work with orthogonal polynomials with degrees well into the millions. From 28273776c3b8d804fe704715432e55971e89f1f9 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 7 May 2024 12:36:14 +0530 Subject: [PATCH 25/25] Remove unused leftover fixkind hack (#240) --- src/chebyshevtransform.jl | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/chebyshevtransform.jl b/src/chebyshevtransform.jl index 7920e6d0..e98f4ed5 100644 --- a/src/chebyshevtransform.jl +++ b/src/chebyshevtransform.jl @@ -22,11 +22,9 @@ ChebyshevTransformPlan{T,kind}(plan::FFTW.r2rFFTWPlan{T,K,inplace,N,R}) where {T # jump through some hoops to make inferrable -_fftKtype(::Val{N}, _...) where N = Vector{Int32} - function plan_chebyshevtransform!(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - ChebyshevTransformPlan{T,1,_fftKtype(Val{N}(), dims...),true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + ChebyshevTransformPlan{T,1,Vector{Int32},true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else ChebyshevTransformPlan{T,1}(FFTW.plan_r2r!(x, FIRSTKIND, dims...; kws...)) end @@ -39,7 +37,7 @@ end function plan_chebyshevtransform(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - ChebyshevTransformPlan{T,1,_fftKtype(Val{N}(), dims...),false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + ChebyshevTransformPlan{T,1,Vector{Int32},false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else ChebyshevTransformPlan{T,1}(FFTW.plan_r2r(x, FIRSTKIND, dims...; kws...)) end @@ -196,7 +194,7 @@ inv(P::IChebyshevTransformPlan{T,1}) where {T} = ChebyshevTransformPlan{T,1}(inv function plan_ichebyshevtransform!(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - IChebyshevTransformPlan{T,1,_fftKtype(Val{N}(), dims...),true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + IChebyshevTransformPlan{T,1,Vector{Int32},true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else IChebyshevTransformPlan{T,1}(FFTW.plan_r2r!(x, IFIRSTKIND, dims...; kws...)) end @@ -208,7 +206,7 @@ end function plan_ichebyshevtransform(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - IChebyshevTransformPlan{T,1,_fftKtype(Val{N}(), dims...),false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + IChebyshevTransformPlan{T,1,Vector{Int32},false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else IChebyshevTransformPlan{T,1}(FFTW.plan_r2r(x, IFIRSTKIND, dims...; kws...)) end @@ -335,7 +333,7 @@ ChebyshevUTransformPlan{T,kind}(plan::FFTW.r2rFFTWPlan{T,K,inplace,N,R}) where { function plan_chebyshevutransform!(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - ChebyshevUTransformPlan{T,1,_fftKtype(Val{N}(), dims...),true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + ChebyshevUTransformPlan{T,1,Vector{Int32},true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else ChebyshevUTransformPlan{T,1}(FFTW.plan_r2r!(x, UFIRSTKIND, dims...; kws...)) end @@ -347,7 +345,7 @@ end function plan_chebyshevutransform(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - ChebyshevUTransformPlan{T,1,_fftKtype(Val{N}(), dims...),false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + ChebyshevUTransformPlan{T,1,Vector{Int32},false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else ChebyshevUTransformPlan{T,1}(FFTW.plan_r2r(x, UFIRSTKIND, dims...; kws...)) end @@ -494,7 +492,7 @@ IChebyshevUTransformPlan{T,kind}(F::FFTW.r2rFFTWPlan{T,K,inplace,N,R}) where {T, function plan_ichebyshevutransform!(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - IChebyshevUTransformPlan{T,1,_fftKtype(Val{N}(), dims...),true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + IChebyshevUTransformPlan{T,1,Vector{Int32},true,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else IChebyshevUTransformPlan{T,1}(FFTW.plan_r2r!(x, IUFIRSTKIND, dims...; kws...)) end @@ -506,7 +504,7 @@ end function plan_ichebyshevutransform(x::AbstractArray{T,N}, ::Val{1}, dims...; kws...) where {T<:fftwNumber,N} if isempty(x) - IChebyshevUTransformPlan{T,1,_fftKtype(Val{N}(), dims...),false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() + IChebyshevUTransformPlan{T,1,Vector{Int32},false,N,isempty(dims) ? NTuple{N,Int} : typeof(dims[1])}() else IChebyshevUTransformPlan{T,1}(FFTW.plan_r2r(x, IUFIRSTKIND, dims...; kws...)) end