<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Stationary GMRF simulation with Discrete Fourier Transformation
Step1:  For benchmarking we will perfom a GF simulation.
Based on non-conditional simulation.
Step2:  Simulation of a temporal GMRF with DFT
Step3:  Now let's build the circulant matrix for the tourus
Oke, for the moment I´ll follow the example in GMRF book.
i.e. a Torus (stationary of 128x128)
Step4:  Algorithm to simulate GMRF with block-circulant
Step5:  Example to perform a FFT in two dimensions
<ASSISTANT_TASK:>
Python Code:
# Load Biospytial modules and etc.
%matplotlib inline
import sys
sys.path.append('/apps')
sys.path.append('..')
#sys.path.append('../../spystats')
import django
django.setup()
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
## Use the ggplot style
plt.style.use('ggplot')
from external_plugins.spystats.spystats import tools as sptools
import scipy
Explanation: Stationary GMRF simulation with Discrete Fourier Transformation
End of explanation
vm = sptools.ExponentialVariogram(sill=0.3,range_a=0.4)
%time xx,yy,z = sptools.simulatedGaussianFieldAsPcolorMesh(vm)
plt.imshow(z)
x = xx[1,:]
y = yy[:,1]
Explanation: For benchmarking we will perfom a GF simulation.
Based on non-conditional simulation.
End of explanation
import scipy.fftpack as fft
c_delta = lambda d : np.hstack(((2 + d),-1,np.zeros(128 - 3),-1))
c_base = c_delta(0.1)
c_base
C = scipy.linalg.circulant(c_base)
n = C.shape[1]
C_inv = np.linalg.inv(C)
plt.imshow(C_inv, interpolation = 'None')
C_inv.max()
plop = fft.ifft(fft.fft(c_base) ** -1)
#plop = fft.rifft(fft.rfft(c_base) ** -1) / n
plop.real.max()
C_inv2 = scipy.linalg.circulant(plop.real)
plt.imshow(C_inv2, interpolation = 'None')
C_chol = np.linalg.cholesky(C_inv)
z = np.random.normal(0, 1, [n])
mmm = np.matmul(C_chol, z)
plt.plot(mmm);
lambda_vec = fft.fft(c_base)
Lambda_aux = np.power(lambda_vec, -0.5)
z_re = np.random.normal(0, 1, [n])
z_im = np.random.normal(0, 1, [n])
z = z_re + 1j * z_im
# z = np.random.normal(0, 1, [n])
x = fft.fft(Lambda_aux.real * z).real / np.sqrt(n)
plt.plot(x);
Explanation: Simulation of a temporal GMRF with DFT
End of explanation
#c_delta = lambda d : np.hstack(((4 + d),-1,np.zeros(128 - 3),-1))
#c_delta = lambda d : np.hstack(((0),-1,np.zeros(128 - 3),-1))
n = 64
N = 64
delta = 0.1
c_base = np.zeros([n, N])
c_base[0, 0] = 4 + delta
c_base[0, 1] = -1
c_base[0, 2] = -1
c_base[1, 0] = -1
c_base[2, 0] = -1
c_base[0, N-1] = -1
c_base[0, N-2] = -1
c_base[N-1, 0] = -1
c_base[N-2, 0] = -1
c_base
%%time 
lambda_mat = fft.fft2(c_base)
z_re = np.random.normal(0, 1, [n, N])
z_im = np.random.normal(0, 1, [n, N])
z = z_re + 1j * z_im
x = fft.fft2((lambda_mat ** -0.5) * z).real / np.sqrt(n *N)
plt.imshow(x, interpolation = 'None')
%%time
n = 64
N = 64
delta = 0.0001
c_base = np.zeros([n, N])
c_base[0, 0] = 4 + delta
c_base[0, 1] = -1
#c_base[0, 2] = -1
c_base[1, 0] = -1
#c_base[2, 0] = -1
c_base[0, N-1] = -1
#c_base[0, N-2] = -1
c_base[N-1, 0] = -1
#c_base[N-2, 0] = -1
 
lambda_mat = fft.fft2(c_base)
z_re = np.random.normal(0, 1, [n, N])
z_im = np.random.normal(0, 1, [n, N])
z = z_re + 1j * z_im
x = fft.fft2((lambda_mat ** -0.5) * z).real / np.sqrt(n *N)
plt.imshow(x, interpolation = 'none')
Explanation: Now let's build the circulant matrix for the tourus
Oke, for the moment I´ll follow the example in GMRF book.
i.e. a Torus (stationary of 128x128)
End of explanation
## Simulate random noise (Normal distributed)
zr = scipy.stats.norm.rvs(size=(C.size,2),loc=0,scale=1)
zr.dtype=np.complex_
#plt.hist(zr.real)
from scipy.fftpack import ifft2, fft2
Lm = scipy.sqrt(C.shape[0]*C.shape[0]) * fft2(C)
Lm.shape
zr.shape
Lm.size
%time v = fft2(scipy.sqrt(Lm) * zr.reshape(Lm.shape))
x = v.real
x.shape
plt.imshow(x,interpolation='None')
cc = scipy.linalg.inv(C)
plt.plot(cc[:,0])
n = x.shape[0]
mm = scipy.stats.multivariate_normal(np.zeros(n),cc)
mmm = mm.rvs()
plt.imshow(mmm.reshape(100,100))
scipy.stats.multivariate_normal?
nn = mm.rvs()
Explanation: Algorithm to simulate GMRF with block-circulant
End of explanation
from scipy.fftpack import ifftn
import matplotlib.pyplot as plt
import matplotlib.cm as cm
N = 30
f, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, sharex='col', sharey='row')
xf = np.zeros((N,N))
xf[0, 5] = 1
xf[0, N-5] = 1
Z = ifftn(xf)
ax1.imshow(xf, cmap=cm.Reds)
ax4.imshow(np.real(Z), cmap=cm.gray)
xf = np.zeros((N, N))
xf[5, 0] = 1
xf[N-5, 0] = 1
Z = ifftn(xf)
ax2.imshow(xf, cmap=cm.Reds)
ax5.imshow(np.real(Z), cmap=cm.gray)
xf = np.zeros((N, N))
xf[5, 10] = 1
xf[N-5, N-10] = 1
Z = ifftn(xf)
ax3.imshow(xf, cmap=cm.Reds)
ax6.imshow(np.real(Z), cmap=cm.gray)
plt.show()
Explanation: Example to perform a FFT in two dimensions
End of explanation
<END_TASK>