Skip to content

Support 3D tensors #185

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Aug 12, 2022
Merged

Support 3D tensors #185

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8665bc4
Start supporting 3D tensors
dlfivefifty Aug 12, 2022
d06e5f0
Tensor Chebyshev 2
dlfivefifty Aug 12, 2022
54919b1
ichebyshevtransform
dlfivefifty Aug 12, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "FastTransforms"
uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
version = "0.14.3"
version = "0.14.4"

[deps]
AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
Expand Down
200 changes: 125 additions & 75 deletions src/chebyshevtransform.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,23 +54,100 @@ plan_chebyshevtransform(x::AbstractArray, dims...; kws...) = plan_chebyshevtrans
@inline _plan_mul!(y::AbstractArray{T}, P::Plan{T}, x::StridedArray{T}) where T = mul!(y, P, x)
@inline _plan_mul!(y::AbstractArray{T}, P::Plan{T}, x::AbstractArray) where T = mul!(y, P, convert(Array{T}, x))

@inline _cheb1_rescale!(_, y::AbstractVector) = (y[1] /= 2; ldiv!(length(y), y))

@inline function _cheb1_rescale!(d::Number, y::AbstractMatrix{T}) where T

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)
ldiv!(2, view(y,1,:))
lmul!(α, @view(y[1,:]))
else
ldiv!(2, view(y,:,1))
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]))
end
end


@inline function _cheb1_rescale!(d::Number, y::AbstractArray)
ldiv_dim_begin!(2, d, y)
ldiv!(size(y,d), y)
end

# TODO: higher dimensional arrays
@inline function _cheb1_rescale!(d::UnitRange, y::AbstractMatrix{T}) where T
@assert d == 1:2
ldiv!(2, view(y,1,:))
ldiv!(2, view(y,:,1))
ldiv!(prod(size(y)), y)
function _prod_size(sz, d)
ret = 1
for k in d
ret *= sz[k]
end
ret
end


@inline function _cheb1_rescale!(d::UnitRange, y::AbstractArray)
for k in d
ldiv_dim_begin!(2, k, y)
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}
Expand All @@ -90,27 +167,21 @@ function mul!(y::AbstractArray{T,N}, P::ChebyshevTransformPlan{T,1,K,false,N}, x
end


_cheb2_rescale!(_, y::AbstractVector) = (y[1] /= 2; y[end] /= 2; ldiv!(length(y)-1, y))

function _cheb2_rescale!(d::Number, y::AbstractMatrix{T}) where T
if isone(d)
ldiv!(2, @view(y[1,:]))
ldiv!(2, @view(y[end,:]))
else
ldiv!(2, @view(y[:,1]))
ldiv!(2, @view(y[:,end]))
end
function _cheb2_rescale!(d::Number, y::AbstractArray)
ldiv_dim_begin!(2, d, y)
ldiv_dim_end!(2, d, y)
ldiv!(size(y,d)-1, y)
end

# TODO: higher dimensional arrays
function _cheb2_rescale!(d::UnitRange, y::AbstractMatrix{T}) where T
@assert d == 1:2
ldiv!(2, @view(y[1,:]))
ldiv!(2, @view(y[end,:]))
ldiv!(2, @view(y[:,1]))
ldiv!(2, @view(y[:,end]))
ldiv!(prod(size(y) .- 1), y)
function _cheb2_rescale!(d::UnitRange, y::AbstractArray)
for k in d
ldiv_dim_begin!(2, k, y)
ldiv_dim_end!(2, k, y)
end

ldiv!(_prod_size(size(y) .- 1, d), y)
end

function *(P::ChebyshevTransformPlan{T,2,K,true,N}, x::AbstractArray{T,N}) where {T,K,N}
Expand Down Expand Up @@ -200,33 +271,25 @@ 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 _icheb1_prescale!(_, x::AbstractVector) = (x[1] *= 2)
@inline function _icheb1_prescale!(d::Number, x::AbstractMatrix)
if isone(d)
lmul!(2, view(x,1,:))
else
lmul!(2, view(x,:,1))
end
@inline function _icheb1_prescale!(d::Number, x::AbstractArray)
lmul_dim_begin!(2, d, x)
x
end
@inline function _icheb1_prescale!(d::UnitRange, x::AbstractMatrix)
lmul!(2, view(x,:,1))
lmul!(2, view(x,1,:))
@inline function _icheb1_prescale!(d::UnitRange, x::AbstractArray)
for k in d
_icheb1_prescale!(k, x)
end
x
end
@inline _icheb1_postscale!(_, x::AbstractVector) = (x[1] /= 2)
@inline function _icheb1_postscale!(d::Number, x::AbstractMatrix)
if isone(d)
ldiv!(2, view(x,1,:))
else
ldiv!(2, view(x,:,1))
end
@inline function _icheb1_postscale!(d::Number, x::AbstractArray)
ldiv_dim_begin!(2, d, x)
x
end

@inline function _icheb1_postscale!(d::UnitRange, x::AbstractMatrix)
ldiv!(2, view(x,1,:))
ldiv!(2, view(x,:,1))
@inline function _icheb1_postscale!(d::UnitRange, x::AbstractArray)
for k in d
_icheb1_postscale!(k, x)
end
x
end

Expand All @@ -249,40 +312,27 @@ function mul!(y::AbstractArray{T,N}, P::IChebyshevTransformPlan{T,1,K,false,N},
ldiv!(2^length(P.plan.region), y)
end

@inline _icheb2_prescale!(_, x::AbstractVector) = (x[1] *= 2; x[end] *= 2)
@inline function _icheb2_prescale!(d::Number, x::AbstractMatrix)
if isone(d)
lmul!(2, @view(x[1,:]))
lmul!(2, @view(x[end,:]))
else
lmul!(2, @view(x[:,1]))
lmul!(2, @view(x[:,end]))
end
@inline function _icheb2_prescale!(d::Number, x::AbstractArray)
lmul_dim_begin!(2, d, x)
lmul_dim_end!(2, d, x)
x
end
@inline function _icheb2_prescale!(d::UnitRange, x::AbstractMatrix)
lmul!(2, @view(x[1,:]))
lmul!(2, @view(x[end,:]))
lmul!(2, @view(x[:,1]))
lmul!(2, @view(x[:,end]))
@inline function _icheb2_prescale!(d::UnitRange, x::AbstractArray)
for k in d
_icheb2_prescale!(k, x)
end
x
end
@inline _icheb2_postrescale!(_, x::AbstractVector) = (x[1] /= 2; x[end] /= 2)
@inline function _icheb2_postrescale!(d::Number, x::AbstractMatrix)
if isone(d)
ldiv!(2, @view(x[1,:]))
ldiv!(2, @view(x[end,:]))
else
ldiv!(2, @view(x[:,1]))
ldiv!(2, @view(x[:,end]))
end

@inline function _icheb2_postrescale!(d::Number, x::AbstractArray)
ldiv_dim_begin!(2, d, x)
ldiv_dim_end!(2, d, x)
x
end
@inline function _icheb2_postrescale!(d::UnitRange, x::AbstractMatrix)
ldiv!(2, @view(x[1,:]))
ldiv!(2, @view(x[end,:]))
ldiv!(2, @view(x[:,1]))
ldiv!(2, @view(x[:,end]))
@inline function _icheb2_postrescale!(d::UnitRange, x::AbstractArray)
for k in d
_icheb2_postrescale!(k, x)
end
x
end
@inline function _icheb2_rescale!(d::Number, y::AbstractArray{T}) where T
Expand All @@ -292,7 +342,7 @@ end
end
@inline function _icheb2_rescale!(d::UnitRange, y::AbstractArray{T}) where T
_icheb2_prescale!(d, y)
lmul!(prod(convert.(T, size(y) .- 1)./2), y)
lmul!(_prod_size(convert.(T, size(y) .- 1)./2, d), y)
y
end

Expand Down
50 changes: 50 additions & 0 deletions test/chebyshevtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -203,6 +203,56 @@ using FastTransforms, Test
@test_throws ArgumentError ichebyshevtransform!(copy(X), Val(2))
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

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))
Expand Down