<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
Step2:  Algorithm to simulate GMRF with block-circulant Matrix.
Taken from
Step3:  For benchmarking we will perfom a GF simulation.
Based on non-conditional simulation.
Step4:  comparison
| Size     | Method        | Seconds   |
|----------|---------------|-----------|
| 100x100  | Full Gaussian | 346       | 
| 100x100  | Markov FFT    | 0.151     |
The stationary circulant markov FFT method is 2291x faster
<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
#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))
#C = scipy.linalg.circulant(c_delta(0.1))
def createToroidalCircularBase(d=0.1,N=128):
    
    Creates a circular base similar to the one described in GMRF Rue and Held, 2005.
    
    c00 = np.hstack(((4 + d),-1,np.zeros(N - 3),-1))
    c01 = np.hstack((-1,np.zeros(N - 1)))
    c0 = np.zeros((N - 2 ,N))
    c1 = np.vstack((c00,c01))
    c = np.vstack((c1,c0))
    c[N -1, 0] = -1
    return c
%%time 
## Create circular base
d = 0.00001
N = 100
c = createToroidalCircularBase(d=d,N=N)
## Simulate random noise (Normal distributed)
from scipy.fftpack import ifft2, fft2
zr = scipy.stats.norm.rvs(size=(c.size,2),loc=0,scale=1,random_state=1234)
zr.dtype=np.complex_
#plt.hist(zr.real)
#Lm = scipy.sqrt(C.shape[0]*C.shape[0]) * fft2(C)
Lm = fft2(c)
v = 1.0/ len(c) * fft2((Lm ** -0.5)* zr.reshape(Lm.shape))
x = v.real
plt.imshow(x,interpolation='None')
## Calculate inverse of c
C_inv = ifft2 ((fft2(c) ** -1))
plt.plot(C_inv[:,0])
Explanation: Algorithm to simulate GMRF with block-circulant Matrix.
Taken from: Rue, H., & Held, L. (2005). Gaussian Markov random fields: theory and applications. CRC press.
Algorithm 2.10
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
%%time 
vm = sptools.ExponentialVariogram(sill=0.3,range_a=0.4)
xx,yy,z = sptools.simulatedGaussianFieldAsPcolorMesh(vm,grid_sizex=100,grid_sizey=100,random_seed=1234)
plt.imshow(z)
Explanation: For benchmarking we will perfom a GF simulation.
Based on non-conditional simulation.
End of explanation
346 / 0.151
c
plt.figure(figsize=(10, 5))
plt.subplot(1,2,1)
plt.imshow(z)
plt.subplot(1,2,2)
plt.imshow(x,interpolation='None')
Explanation: comparison
| Size     | Method        | Seconds   |
|----------|---------------|-----------|
| 100x100  | Full Gaussian | 346       | 
| 100x100  | Markov FFT    | 0.151     |
The stationary circulant markov FFT method is 2291x faster
End of explanation
<END_TASK>