SPDE for Markovian Model¶

In [1]:
import gstlearn as gl
import gstlearn.plot as gp
import numpy as np
import matplotlib.pyplot as plt
In [2]:
Markov = True
In [3]:
dx = [0.5, 0.5]
nx = [400, 400]
db = gl.DbGrid.create(nx, dx)
In [4]:
mesh = gl.MeshETurbo(db, False)
In [5]:
scale = 1
sill = 5.0
param = 2.0

ranges = [1.8 * scale, 1.0 * scale]
angles = [30.0, 0]
if Markov:
    model = gl.Model.createFromParam(
        gl.ECov.MARKOV, ranges=ranges, angles=angles, sill=sill, flagRange=False
    )
    cova = model.getCovAniso(0)
    cova.setMarkovCoeffsBySquaredPolynomials([1, -1.5], [1, -1], 0.001)
    # coeffs = np.array(cova.getMarkovCoeffs())
    # for i in range(len(coeffs)):
    #    coeffs[i]*= scale**(2*i)
    # cova.setMarkovCoeffs(coeffs)
else:
    model = gl.Model.createFromParam(
        gl.ECov.MATERN,
        param=param,
        ranges=ranges,
        sill=sill,
        angles=angles,
        flagRange=False,
    )

model
Out[5]:
Model characteristics
=====================
Space dimension              = 2
Number of variable(s)        = 1
Number of basic structure(s) = 1
Number of drift function(s)  = 0
Number of drift equation(s)  = 0

Covariance Part
---------------
Markov (Third Parameter = 1)
- Sill         =       5.000
- Ranges       =      12.471      6.928
- Theo. Ranges =       1.800      1.000
- Angles       =      30.000      0.000
- Rotation Matrix
                 [,  0]     [,  1]
      [  0,]      0.866     -0.500
      [  1,]      0.500      0.866
Total Sill     =       5.000
Known Mean(s)      0.000
In [6]:
Q = gl.PrecisionOp(mesh, model.getCovAniso(0))
In [7]:
indx = int(nx[0] / 2)
indy = int(nx[1] / 2)
ind = indy + nx[1] * indx
covy = np.array(Q.computeCov(ind)).reshape(nx)
In [8]:
xx = np.arange(len(covy) / 2) * dx[0]
if not Markov:
    covM = [model.evalCov(gl.SpacePoint([0, 0]), gl.SpacePoint([i, 0])) for i in xx]
    plt.plot(xx, covM, "--", label="Covariance")
plt.plot(
    xx, covy[indx, indy : (indy + nx[1])], label="discretized"
)  # /covy[indx,indx])
covy[indx, indx]
ax = plt.legend()
plt.show()
i = plt.imshow(covy.reshape(nx), origin="lower")
No description has been provided for this image
No description has been provided for this image
In [9]:
dbg = gl.DbGrid.create(nx, dx)
x = Q.simulateOne()
dbg["simu"] = x

res = gp.plot(dbg)
aa = plt.axis("equal")
No description has been provided for this image
In [10]:
N = 2**8
cova = model.getCovAniso(0)
mH1 = 3.5 * np.max(cova.getRanges())
mH2 = 3.5 * np.max(cova.getRanges())

result = np.array(cova.evalCovFFT([mH1, mH2], N).getValues())
res = np.array(result).reshape((N, N))

i = plt.imshow(res, extent=[-mH1, mH1, -mH2, mH2], origin="lower")
c = plt.colorbar()
No description has been provided for this image
In [11]:
X1 = np.linspace(-mH1, mH1, N)
X2 = np.linspace(-mH2, mH2, N)
plt.plot(X1, res[int(N / 2), :], label="FFT")

uu = covy[indx, indy : (indy + 100)]
plt.plot(np.arange(len(uu)) * dx[0], uu, "--", label="Discretized covariance")
ax = plt.legend()
No description has been provided for this image
In [12]:
print(round(res[int(N / 2), int(N / 2)] / uu[0], 3))

1.009