I am interested in an efficient code to generate an $n$-D Gaussian random field (sometimes called processes in other fields of research), which has applications in cosmology.
Attempt
I wrote the following code:
fftIndgen[n_] := Flatten[{Range[0., n/2.], -Reverse[Range[1., n/2. - 1]]}]
Clear[GaussianRandomField];
GaussianRandomField::usage = "GaussianRandomField[size,dim,Pk] returns
a Gaussian random field of size size (default 256) and dimensions dim
(default 2) with a powerspectrum Pk";
GaussianRandomField[size_: 256, dim_: 2, Pk_: Function[k, k^-3]] :=
Module[{noise, amplitude, Pk1, Pk2, Pk3, Pk4},
Which[
dim == 1,Pk1[kx_] :=
If[kx == 0 , 0, Sqrt[Abs[Pk[kx]]]]; (*define sqrt powerspectra*)
noise = Fourier[RandomVariate[NormalDistribution[], {size}]]; (*generate white noise*)
amplitude = Map[Pk1, fftIndgen[size], 1]; (*amplitude for all frequels*)
InverseFourier[noise*amplitude], (*convolve and inverse fft*)
dim == 2,
Pk2[kx_, ky_] := If[kx == 0 && ky == 0, 0, Sqrt[Pk[Sqrt[kx^2 + ky^2]]]];
noise = Fourier[RandomVariate[NormalDistribution[], {size, size}]];
amplitude = Map[Pk2 @@ # &, Outer[List, fftIndgen[size], fftIndgen[size]], {2}];
InverseFourier[noise*amplitude],
dim > 2, "Not supported"]
]
Here are a couple of examples on how to use it in one and 2D
GaussianRandomField[1024, 1, #^(-1) &] // ListLinePlot
GaussianRandomField[] // GaussianFilter[#, 20] & // MatrixPlot
Question
The performance is not optimal — On other interpreted softwares, such Gaussian random fields can be generated ~20 times faster. Do you have ideas on how to speed things up/improve this code?
AbsoluteTiming
andPrint
statements tell you? $\endgroup$size = 1024*1024; Pk1[k_] := If[k == 0 , 0, Sqrt[Abs[k^(-1)]]]; (*define sqrt powerspectra*);
noise = Fourier[ RandomVariate[ NormalDistribution[], {size}]]; (*generate white noise*)// AbsoluteTiming
amplitude = Map[Pk1, fftIndgen[size], 1]; (*amplitude for all frequels*)// AbsoluteTiming
InverseFourier[noise*amplitude]; // AbsoluteTiming
So I guess the Map Powerpsectra part ? {0.194598,Null} {4.388953,Null} {0.506372,Null} $\endgroup$