From 2dcabebd3b31a59b3ece09e390e127a279a941c8 Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Sat, 21 Oct 2023 21:41:30 +0100 Subject: [PATCH] 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