Power spectrum for n-dimensional LowerTriangularArray

[milankl]

This allows to calculate the power spectrum on n-dimensional LowerTriangularArrays along the first/spherical harmonic dimension, e.g.

julia> L
15×2 LowerTriangularArray{ComplexF32, 2, Matrix{ComplexF32}}:
 -0.382899-0.119496im   -0.403847+0.663406im
 -0.364989+0.0916957im   0.258135+0.630038im
 -0.712897-0.499565im    -1.26181-0.6088im
 -0.700691+0.644983im   -0.360479+0.602015im
   1.12504+0.768878im   -0.579596+1.31353im
  -1.29236-0.566898im    0.688637+0.0492758im
 -0.732031-1.05859im      1.06551-1.45981im
 -0.601207+0.140327im     0.61631+0.334235im
  0.244479-0.205514im     0.38642+0.428516im
 0.0773385+1.09861im     -1.01195-1.18484im
  0.950894-0.682852im   -0.301167-0.599011im
 -0.628636+0.487482im    0.450988+0.312973im
 -0.331025+1.8114im     -0.648268+0.942647im
   1.26967+0.858582im   -0.424846-0.124433im
  0.107177-0.159691im   -0.127849-0.105034im

julia> L[:, 1]
15-element, 5x5 LowerTriangularMatrix{ComplexF32}
 -0.382899-0.119496im         0.0+0.0im       …        0.0+0.0im            0.0+0.0im
 -0.364989+0.0916957im   -1.29236-0.566898im           0.0+0.0im            0.0+0.0im
 -0.712897-0.499565im   -0.732031-1.05859im            0.0+0.0im            0.0+0.0im
 -0.700691+0.644983im   -0.601207+0.140327im     -0.331025+1.8114im         0.0+0.0im
   1.12504+0.768878im    0.244479-0.205514im       1.26967+0.858582im  0.107177-0.159691im

julia> SpeedyTransforms.power_spectrum(L)
5×2 Matrix{Float32}:
 0.160891  0.6032
 1.37493   0.472293
 1.29932   2.67026
 1.59882   0.713167
 0.899884  0.41962

julia> SpeedyTransforms.power_spectrum(L[:, 1])
5-element Vector{Float32}:
 0.16089073
 1.374925
 1.2993188
 1.5988191
 0.8998841

julia> SpeedyTransforms.power_spectrum(L[:, 2])
5-element Vector{Float32}:
 0.6032003
 0.4722927
 2.6702583
 0.7131671
 0.4196203

[milankl]

@maximilian-gelbrecht I've added this method because L[l, m, k] is quite a common access pattern arising in

    for k in eachmatrix(spec)
        for m in 1:mmax
            for l in m:lmax 
                L[l, m, k]

to loop over integers l, m but any additional dimensions k (eachmatrix loops over all of those, returning a CartesianIndex. I added the type contraints of i::Integer, j::Integer, c::CartesianIndex to avoid any dispatch issues, but not sure they're actually needed, e.g. i, j, c::CartesianIndex could work equally well.

[maximilian-gelbrecht]

I think i, j, c::CartesianIndex might cause some problems with the getindex functions that use Vararg. I could be mistaking though.

Anyway, good to a function like this

[milankl]

Turns out this method wasn't needed, but we had a @boundscheck M == N+1 in getindex with varargs which you added to guarantee that someones wants to do matrix indexing with e.g. 2 indices on a LowerTriangularMatrix. But checking this from the number of indices isn't helpful because L[1, 1, 1, 1, 1, 1] is the same as L[1] regardless the shape of L and also you can add many trailing CartesianIndex() that would increase the number of indices even though they don't actually count. Just removed that @boundscheck to see what happens, I've tested it manually it doesn't seem necessary!