I need to prepare "flattened" versions of 2D fftfrequencies in the shape Nx^2 * 2. Those are basically constructed like a ravel(meshgrid(fftfreqs1d,fftfreqs1d)) in matlab or python.

This appears to be no big deal in python, but can hang for reasonable array sizes in julia, especially when i want to build a StaticArray out of the intermediate results. To make it more confusing, @btime pretends that my arrays are created in no time, while they are clearly not. My question is why this happens and how it is done right.

I am aware that using julia it might be a waste to keep the full 2D fftfreqs in memory instead of using the 1D versions and a loop, but let us assume for a moment that i need it this way.

Julia

function my_freqs1(Nnu::Int,T)
    dx = 2. /Nnu
    freq1d = fftfreq(Nnu).*dx
    nu = hcat(  vec([ i for i in freq1d, j in freq1d  ]),
                                    vec([ j for i in freq1d, j in freq1d  ]))
    return nu
end;
@btime my_freqs1(100,Float64)

28.528 μs (10 allocations: 312.80 KiB)

Julia, converting to a static array (in the hope for better performance of other code later on)

function my_freqs2(Nnu::Int,T)
    ### the same as above ###
    return SMatrix{Nnu^2,2,T}(nu)
end;
@btime my_freqs2(100,Float64)

94.540 μs (36 allocations: 470.38 KiB)

Python

def my_fftfreqs(xy):
    freqs = np.fft.fftfreq(np.shape(xy)[0],d=xy[1]-xy[0])
    fx,fy = np.meshgrid(freqs,freqs,indexing="ij")
    freq_list = np.transpose(np.asarray( [np.ravel(fx),np.ravel(fy)] ))
    return freq_list
%time f=my_fftfreqs(np.linspace(0,1,100));

CPU times: user 1.08 ms, sys: 0 ns, total: 1.08 ms
Wall time: 600 µs

My observation is that while python %time reports a much longer time, it will actually run in a very reasonable time while the julia version has a noticable delay and the version with the static array will hang for a long time and completely crash for larger sizes.

Please help me to understand how i would do this correctly in Julia and whether (why not?) creating a static array seems to be such a bad idea.

CodePudding user response：

Rather than making a SMatrix{Nnu^2,2} I think you probably want to make a Vector{SVector{2}}. The former will require recompiling for each new value of Nnu which is fairly inefficient.

CodePudding user response：

You may also consider:

using FFTW

my_freqs3(ν) = fftfreq(ν)*2/ν |> 
  (w -> [repeat(w, inner=length(w)) repeat(w, outer=length(w))])

# or

my_freqs3alt(ν) = ( w = fftfreq(ν)*2/ν ; 
  [repeat(w, inner=length(w)) repeat(w, outer=length(w))] )

which is more Julian and "if-I-understand-correctly" is equivalent.

Usually shorter/simpler functions are also more efficient.

Julia features used:

• Unicode nu variable.
• Piping |> operator.
• Definition with no function keyword.
• repeat standard library vector filling function.
• Matlab-like hcat [v1 v2] notation.
• Multi-statement block enclosed in ( ) separated by ;.