Skip to content
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.

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

Type inconsistencies with FFTW? #14

Open
milankl opened this issue Sep 1, 2022 · 6 comments
Open

Type inconsistencies with FFTW? #14

milankl opened this issue Sep 1, 2022 · 6 comments

Comments

@milankl
Copy link
Collaborator

milankl commented Sep 1, 2022

FFTW plans are <:AbstractFFTs.Plan{T} with T complex for complex-to-real plans and T real for real-to-complex plans. So the parametric type of the plan is the type of the input. However, the plans here always seem to be <:AbstractFFTs.Plan{Complex}. Is that on purpose?

julia> generic_plan = GenericFFT.plan_rfft(zeros(Float16,64))
GenericFFT.DummyrFFTPlan{ComplexF16, false, UnitRange{Int64}}(64, 1:1, #undef)

julia> fftw_plan = FFTW.plan_rfft(zeros(64))
FFTW real-to-complex plan for 64-element array of Float64
(rdft2-ct-dit/8
  (hc2c-direct-8/28/1 "hc2cfdftv_8_avx2"
    (rdft2-r2hc-direct-8 "r2cf_8")
    (rdft2-r2hc01-direct-8 "r2cfII_8"))
  (dft-direct-8-x4 "n2fv_8_avx2_128"))

julia> supertype(supertype(typeof(fftw_plan)))
AbstractFFTs.Plan{Float64}

julia> supertype(supertype(typeof(generic_plan)))
AbstractFFTs.Plan{ComplexF16}     # expected Float16 not ComplexF16
@milankl
Copy link
Collaborator Author

milankl commented Sep 1, 2022

I believe these lines

GenericFFT.jl/src/fft.jl

Lines 305 to 306 in 926042f

plan_rfft(x::StridedArray{T}, region) where {T <: RealFloats} = DummyrFFTPlan{Complex{real(T)},false,typeof(region)}(length(x), region)
plan_brfft(x::StridedArray{T}, n::Integer, region) where {T <: ComplexFloats} = DummybrFFTPlan{Complex{real(T)},false,typeof(region)}(n, region)

could be changed to

plan_rfft(x::StridedArray{T}, region) where {T <: RealFloats} = DummyrFFTPlan{T,false,typeof(region)}(length(x), region) 
plan_brfft(x::StridedArray{T}, n::Integer, region) where {T <: ComplexFloats} = DummybrFFTPlan{T,false,typeof(region)}(n, region) 

as where already enforces real/complex T. Happy to open a PR just wanted to check that's fine?

@daanhb
Copy link
Member

daanhb commented Sep 1, 2022

I don't know about the history here (perhaps @MikaelSlevinsky does) but it seems better to me to mimick what FFTW does. Thanks for catching this. Yes, we'll definitely consider a PR.

This would be a breaking change I guess.

@MikaelSlevinsky
Copy link
Member

Not sure about the dummy plans. Maybe @dlfivefifty knows? Probably it's best to do as you say.

@daanhb
Copy link
Member

daanhb commented Sep 1, 2022

If we're going to change the plans in a (slightly) breaking way, I would like to go for including the size of the arrays in the plans also. That makes the plans less dummy, but it would (a) prevent errors and (b) allow progress with multidimensional transforms.

But one step at a time is also good.

@milankl
Copy link
Collaborator Author

milankl commented Sep 2, 2022

By "size" you mean actual size and not just length? Because the length is already included no?

julia> import GenericFFT

julia> rfft_plan = GenericFFT.plan_rfft(zeros(Float16,1024))
GenericFFT.DummyrFFTPlan{ComplexF16, false, UnitRange{Int64}}(1024, 1:1, #undef)

julia> rfft_plan.n
1024

julia> brfft_plan = GenericFFT.plan_brfft(zeros(ComplexF16,513),1024)
GenericFFT.DummybrFFTPlan{ComplexF16, false, UnitRange{Int64}}(1024, 1:1, #undef)

julia> brfft_plan.n
1024

@daanhb
Copy link
Member

daanhb commented Sep 2, 2022

Good point. It seems the length is included in some plans, but not in all. It is included for these plans here with an Integer field.

I had in mind fixing issue #10. With FFTW it is like this:

julia> using FFTW

julia> p = FFTW.plan_rfft(zeros(1024));

julia> p.sz
(1024,)

julia> p = FFTW.plan_rfft(zeros(16,16));

julia> p.sz
(16, 16)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants