# This source code is public domain 
# Author: Christian Schirm
import numpy
import matplotlib.pyplot as plt

# To speed up the random generation, a numerical trick 
# is used based on the fast fourier transform. See e.g.
# D. P. Kroese et al., "Spatial Process Generation", 2013
# https://arxiv.org/pdf/1308.0399.pdf

n = 2**8
x = numpy.arange(-n/2,n/2)/float(n)
def covFunc(r): return numpy.exp(-r)
def cov(deltaX, deltaY): return covFunc(8*(deltaX**2 + deltaY**2)**0.5)
indX,indY = numpy.mgrid[0:n,0:n]
C = cov(x[indX],x[indY])
Gamma = numpy.fft.fft2(numpy.roll(numpy.roll(C,len(C)/2,0),len(C)/2,1))
Z1 = numpy.random.randn(n,n) + 1J*numpy.random.randn(n,n)
X = numpy.fft.fft2(numpy.sqrt(Gamma) *Z1/n).real

fig, (plt2, plt1) = plt.subplots(nrows=1, ncols=2, figsize=(10./1.7, 4.5/1.7), dpi=120)
plt1.imshow(X,cmap='viridis',vmin=-2.5,vmax=2.5,extent=[0,10,0,10])
xCov = numpy.linspace(0,0.825,40)
plt2.plot(8*xCov, covFunc(8*xCov), label='exp(-|r|)')
plt2.set_xlabel('r')
plt2.set_ylabel('Covariance')
plt2.axis([0,5,0,1])
plt2.legend()
plt.savefig('Gaussian_process_2D_abs_exp.png')