Autovectorization and AbstractFFTs

[andreasnoack]

I'm opening this issue to discuss something I started thinking about when working on #76. I have mixed feelings about ffts of lower dimension than the array. It feels like "autovectorization" which is generally avoided in other parts of Julia's multidimensional code. It would be more in line with the general array code to apply a lower-dimensional fft plan along some dimensions with something like mapslices.

Maybe some of this is imposed by AbstractFFTs, but I'm also now in doubt what this package really gains from using that interface. Some of the promotion methods might be convenient, but it isn't that big of a deal IMO, and it essentially forces the interface of FFTW while also introducing issues like #66. AbstractFFTs doesn't solve the problem of switching backends for downstream packages (although there is ongoing work that might change that), so I'm not sure if depending on AbstractFFTs is a net gain. My main goal with working on this package is to allow e.g. https://github.com/JuliaStats/KernelDensity.jl to work without FFTW, and there might be easier ways to allow for that than making changes to AbstractFFTs.

[dannys4]

FWIW I don't really like AbstractFFTs for the exact reason you say of "making everything FFTW". If you're not extremely careful, it reallocates your input into a new array which may be significantly less performant---I forget if I put this in an issue or what, but at some point I mentioned that it will take a SVec and just copy it to a new Vector and run the FFT on that. Really irritating to me.

[dannys4]

This is irritating if you're working with SVec, but "could be" debilitating if you're working with, e.g., distributed data. That being said, of course, an often ignored mantra of a performance engineer is "optimize for the use-case you have, not the use-case you could have".

Anyway, all this to say that I'm supportive of whatever decision you make down the line.

[andreasnoack]

For those who might not know what I'm proposing to the automatic vectorization that FFTW defines, here is an example of what I'd consider more idiomatic Julia using mapslices

julia> x = randn(4,5)
4×5 Matrix{Float64}:
  1.71268   -0.271063    -2.09614  0.114768  -0.552492
  0.588278   1.18662     -2.07875  0.840141   0.562409
 -0.862501  -0.625797     1.4668   0.905327   0.3185
 -0.771088   0.00501239  -0.3175   1.34392   -0.406535

julia> p1 = plan_fft(x[:, 1]);

julia> p2 = plan_fft(x[1, :]);

julia> mapslices(fft, x; dims = 1) == mapslices(Base.Fix1(*, p1), x; dims = 1) == fft(x, 1)
true

julia> mapslices(fft, x; dims = 2) == mapslices(Base.Fix1(*, p2), x; dims = 2) == fft(x, 2)
true

[jusack]

I get that the implicit vectorization of FFTW might not be the most Julian way, but I appreciate that it is rather simple to use.
For example

A = randn(5,5)
B = mapslices(rfft, A, dims=1)

A2 = mapslices(Base.Fix2(irfft, 5), B, dims=1)

is a lot longer than

A2 = irfft(B, 5, 1)

or is there a better way?

Edit:
Maybe that idea is a bit crazy, but if the AbstractFFTs interface is dropped and the FFTA interface is simplified to remove all autovectorization. The above could be simplified by adding a one argument version of irfft that is just this call to Base.Fix2, i.e. FFTA.irfft(N::Int) = Base.Fix2(irfft,N) Then the above would write as mapslices(irfft(5), B, dims=1). N could also be a keyword argument

[dannys4]

As it stands, there is some dispatching needed (possibly piracy?) because vanilla Julia B = mapslices(rfft, A, dims=1) will create a callgraph for every slice (Correct me if I'm wrong though). What you would actually want is something that looks like B = mapslices(Base.Fix1(*, FFTA.plan_rfft(A,dims=1)), A, dims=1) (this doesn't work, don't try it).

Anyway, if we're creating a dispatch of mapslices anyway, we could add a new positional or keyword argument....

[jusack]

I just tested it and you are right: mapslices with (r)fft will create a new callgraph for every slice... If the AbstractFFTs interface is dropped I would suggest to have dims as a keyword argument, similar to functions like mean or sum that can also directly be applied along a specific dimension to have efficient vectorization support built-in to the high-level methods.