Extend FTPlans to Arrays and higher dimensions

[ioannisPApapadopoulos]

No description provided.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 83.07%. Comparing base (b2cd0e8) to head (16b6f43).
Report is 15 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #252      +/-   ##
==========================================
+ Coverage   82.78%   83.07%   +0.29%     
==========================================
  Files          16       17       +1     
  Lines        2805     2930     +125     
==========================================
+ Hits         2322     2434     +112     
- Misses        483      496      +13     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[ioannisPApapadopoulos]

@dlfivefifty I've implemented the n-dimensional array * and \ for FTPlan as discussed in JuliaApproximation/ClassicalOrthogonalPolynomials.jl#206 (comment)

Works much faster than looping over columns naively.

[dlfivefifty]

Cool! Note I implemented something very similar (but assuming its a matrix, not a 1D plan) here:

https://github.com/JuliaApproximation/ContinuumArrays.jl/blob/e178b303c8355d88fbb833e578c8c7d74ad75948/src/plans.jl#L45

Probably we should consolidate these two types?

[ioannisPApapadopoulos]

Ah yep, I was not aware. Am I right in thinking MulPlan(matrix, dims) loops over the dimensions via the for loops? If that is the case then that really slows things down with an FTPlan. e.g.

c = randn(5,20,20,20);
F = plan_cheb2leg(c)
P = NDimsPlan(F, size(c), (1,))

@time for i in axes(c,2), j in axes(c,3), k in axes(c,4)
    F*view(c,:,i,j,k)
end # 0.003361 seconds (32.42 k allocations: 1.715 MiB)
# 10x faster via NDimsPlan
@time P*c; # 0.000542 seconds (42 allocations: 1.833 MiB)

[ioannisPApapadopoulos]

Or is it only looping over the dims..?

[jishnub]

Why the array construction in the indexing? Shouldn't this just be P.szs[k...]?

[ioannisPApapadopoulos]

I've probably misunderstood what k is. I copied the size(P, k...) structure from size(P::JacobiTransformPlan, k...) = size(P.chebtransform, k...) and assumed it is a tuple. In which case P.szs[k...] does not work.

e.g.

julia> szs = (10,20,30);
julia> szs[[(1,2)...]]
(10, 20)
julia> szs[(1,2)...]
ERROR: MethodError: no method matching getindex(::Tuple{Int64, Int64, Int64}, ::Int64, ::Int64)

[jishnub]

Perhaps I've misunderstood this. Are you trying to obtain multiple indices in one go?

[ioannisPApapadopoulos]

That's the intention! Although I don't need that functionality in this PR.

[jishnub]

Perhaps best to remove it then, as it's not the conventional usage of size. One may broadcast size over indices to obtain the sizes along multiple dimensions at once.

[dlfivefifty]

Note sizes of plans and arrays behave differently

compare w FFT

[ioannisPApapadopoulos]

judging by this line here https://github.com/JuliaMath/FFTW.jl/blob/f888022d7a1ff78491abf8f33f1055cc52a68f0a/src/fft.jl#L254 I think I should only really keep size(P::NDimsPlan) = P.szs

[ioannisPApapadopoulos]

...or maybe I should keep size(P::NDimsPlan, k::Int) = P.szs[k]?

[ioannisPApapadopoulos]

Any opposition to merging? @dlfivefifty @jishnub

[ioannisPApapadopoulos]

Hi @dlfivefifty @jishnub @MikaelSlevinsky unless someone is opposed, I will bump the version number, merge and tag in a couple of days!

[jishnub]

I don't have any objections

[MikaelSlevinsky]

It does what you want it to do?

[ioannisPApapadopoulos]

Yes! It makes it easy to take an FTPlan and turn it into an n-dimensional FTPlan applied across multiple array dimensions.