Add a Fourier generator for periodic spatial random fields

[LSchueler]

This is a first attempt to include a new generator based on the Fourier method. One of its desirable properties is the periodicity of the generated fields, which can be a necessity for lower boundary conditions in atmospherical physics.

Tests are still completely missing and I think the length scales of the fields are currently not working as expected, which is probably due to the rescaling of the Fourier modes. Before commiting any more to this, I'll wait for the comments from an atmospherical physicist.

Here is a simple script to show you how to currently use the new generator:

import numpy as np
import gstools as gs
import matplotlib.pyplot as pt

x = np.arange(100)
y = np.arange(100)

model = gs.Gaussian(dim=2, var=1, len_scale=30.0)
srf = gs.SRF(model, generator='Fourier', modes_truncation=[1.0, 1.0], modes_no=[80, 80], seed=3684903)

field = srf((x, y), mesh_type='structured')

srf.plot()
pt.show()

[mschlutow]

Something's not quite right. The shape of the histogram is not Gaussian and its variance not unity. Becomes even more crooked when more modes are added.

import numpy as np
import gstools as gs
import matplotlib.pyplot as pt

Nx=512
Ny=256

x = np.arange(Nx)
y = np.arange(Ny)

model = gs.Gaussian(dim=2, var=1, len_scale=100.0)
srf = gs.SRF(model, generator='Fourier', modes_truncation=[1.0, 1.0], modes_no=[80, 80], seed=3684903)

field = srf((x, y), mesh_type='structured')

srf.plot()
pt.show()

pt.hist(field.flatten('C'),bins=64)
pt.show()

[LSchueler]

Thanks for providing the script 👍 I guess the reason is that you only sample very few correlation lengths. I just ran your script, but with len_scale=1.0 and this is the result:

[mschlutow]

Thanks for the quick reply.
I am afraid that it is not only the length scale that leads to uneven distribution.

model = gs.Gaussian(dim=2, var=1.0, len_scale=200.0)
srf = gs.SRF(model, generator='Fourier', modes_truncation=[1.0, 1.0], modes_no=[512, 256])

Why is it symmetrical?

[LSchueler]

Yeah, there is definitely something wrong with the length scales. Here is a script to estimate the variograms, which also gives you the actual length scale of a field. If you want to use your own field generator, simply ignore everything until line 21 and give your field to the variogram estimator as the field keyword in line 29.
I have to see, if I find some time this evening to check the rescaling of the Fourier modes.

import numpy as np
import gstools as gs
import matplotlib.pyplot as pt

x = np.arange(0, 100, 1)
y = np.arange(0, 60, 1)

dim = 2
model = gs.Gaussian(dim=dim, var=1, len_scale=[5.0, 3.0])
# fourier method
srf = gs.SRF(
    model,
    generator='Fourier',
    modes_truncation=[1.0, 1.0],
    modes_no=[100, 60],
    seed=3684903,
)
# randomization method
# srf = gs.SRF(model, seed=3684903)

field = srf((x, y), mesh_type='structured')

x_max = 40.0
angle = 0.0

bins = np.arange(0, x_max, 2)

bin_center, dir_vario, counts = gs.vario_estimate(
    *((x, y), field, bins),
    direction=gs.rotated_main_axes(dim, angle),
    angles_tol=np.pi / 16,
    bandwidth=8,
    mesh_type="structured",
    return_counts=True,
)

model.fit_variogram(bin_center, dir_vario)
print(f'Fitted model: {model}')

fig, (ax1, ax2) = pt.subplots(1, 2, figsize=[10, 5])

ax1.scatter(bin_center, dir_vario[0], label="emp. vario in x-dir")
ax1.scatter(bin_center, dir_vario[1], label="emp. vario: in y-dir")
ax1.legend(loc="lower right")

model.plot("vario_axis", axis=0, ax=ax1, x_max=x_max, label="fit in x-dir")
model.plot("vario_axis", axis=1, ax=ax1, x_max=x_max, label="fit in y-dir")
ax1.set_title("Fitting an anisotropic model")

srf.plot(ax=ax2)
pt.show()

[LSchueler]

Yeah, I think the final problem to solve is the scaling of the Fourier modes. This is an example, where I manually figured out the rescaling for one specific set of parameters and the (anisotropic) length scales are estimated correctly.

[MuellerSeb]

Please rebase with main.

[MuellerSeb]

Thanks @LSchueler for the great work. Will have a look as soon as I find the time!
Could you update the branch with main?

[MuellerSeb]

I didn't test it yet and just had a first look. There are some minor things already but it looks good overall! Great work.

[LSchueler]

I merged the main branch and I also implemented the Fourier method in GSTools-Core, so included that import in GSTools, if GSTools-Core is available.

[MuellerSeb]

This should be solved with #343

[LSchueler]

I'll solve that in a different PR.

[LSchueler]

That's going to be solved in PR #354.

[LSchueler]

PR #354 should be merged first and then I can update this branch with the fixes.

[MuellerSeb]

I tested it and think this is fine! Some comments about the Generator.

One thing I noticed is that variance is too low when estimated but I think this is fine, since ergodicity is hard/impossible to achieve with periodic fields and that it is in the nature of things that field variance is lower than the defined variance of the used model.

My observation is: the lower the period, the lower the variance. And this makes perfectly sense.

[LSchueler]

This is ready for the next review round.

[MuellerSeb]

Like it! some minor comments.

[MuellerSeb]

LGTM

[LSchueler]

I can't believe it, after nearly 1 1/2 years, this is ready! Thanks for the reviews.