SPDE¶

In this tutorial, we show how to use the API SPDE.

In [1]:
import gstlearn as gl
import gstlearn.plot as gp
import numpy as np
import matplotlib.pyplot as plt

from sksparse.cholmod import cholesky
import scipy as sc
from scipy.sparse import *
from scipy.sparse.linalg import *
import numpy as np
In [2]:
# Data
np.random.seed(123)
ndat = 1000

# Model

rangev = 0.2
sill = 1.
nugget = 0.1

# Grid 
nx = [50,50]
dx = [0.02,0.02]
x0 = [0,0]

#Grid meshing

nxm = [75,75]
dxm = [0.02,0.02]
x0m = [-0.25,-0.25]

Grids definition¶

In [3]:
grid = gl.DbGrid.create(nx,dx,x0)
gridExt = gl.DbGrid.create(nxm,dxm,x0m)
mesh = gl.MeshETurbo(gridExt)

Model definition¶

In [4]:
model = gl.Model.createFromParam(gl.ECov.BESSEL_K,param=1,range=rangev,sill=sill)
model.addCovFromParam(gl.ECov.NUGGET,sill=nugget)

SPDE simulation¶

In [5]:
spde = gl.SPDE(model,grid,None,gl.ESPDECalcMode.SIMUNONCOND)

Grid query¶

In [6]:
iuid = spde.compute(grid)
In [7]:
ax = grid.plot()
ax.decoration(title="Non Conditional Simulation")

Data query¶

In [8]:
dat = gl.Db.create()
dat["x"] = np.random.uniform(size=ndat)
dat["y"] = np.random.uniform(size=ndat)
dat.setLocators(["x","y"],gl.ELoc.X)
iuid = spde.compute(dat)

Kriging¶

In [9]:
spdeRes = gl.SPDE(model,grid,dat,gl.ESPDECalcMode.KRIGING,mesh)
In [10]:
iuid = spdeRes.compute(grid)
ax = grid.plot()
ax.decoration(title="Estimation")

Manually¶

Aproj mesh to grid¶

In [11]:
Pglg = gl.ProjMatrix(grid,mesh)
Aprojg = Pglg.getAprojToTriplet().toTL()

Simulation¶

In [12]:
Q = spdeRes.getPrecisionOp().getQToTriplet().toTL()
cholQ =  cholesky(Q)
In [13]:
u = np.random.normal(size = Q.shape[0])
gridExt["simuManual"] = cholQ.apply_Pt(cholQ.solve_Lt(1./np.sqrt(cholQ.D())*u))
gridExt.addSelection((gridExt["x1"]>0) & (gridExt["x2"]>0) & (gridExt["x1"]<1.) & (gridExt["x2"]<1.))
ax = gridExt.plot("simuManual",useSel=False)
round(np.var(gridExt["simuManual"][np.where(gridExt["NewSel"]==1)]),4)
Out[13]:
1.0167

Kriging¶

In [14]:
Pgl = gl.ProjMatrix(dat,mesh)
Aproj = Pgl.getAprojToTriplet().toTL()

Qop = Q + 1/nugget * Aproj.T @ Aproj
cholQop =  cholesky(Qop)

kriging = cholQop.solve_A(Aproj.T @ (dat["spde*"]/nugget))
In [15]:
grid["manually"] = Aprojg @ kriging
In [16]:
ax = plt.scatter(grid["manually"],grid["*kriging"],s=1)
p = plt.plot([-3,3],[-3,3],c="r")

Likelihood¶

Manually with Cholesky vs. matrix-free approach with SPDE api.

In [17]:
def solveMat(cholQop,x):
    return cholQop.solve_A(x)

def invSigma(sigma2,Aproj,cholQop,x):
    return 1./sigma2 * (x - 1./sigma2 * Aproj @ solveMat(cholQop, Aproj.T @ x))

def detQ(cholQ):
    return cholQ.logdet()

x = dat["spde.simu"]
ones = np.ones_like(x)
invSigmaOnes = invSigma(nugget,Aproj,cholQop,ones)
mu  = np.sum(x * invSigmaOnes)/np.sum( ones * invSigmaOnes) 
quad = np.sum((x-mu)*invSigma(nugget,Aproj,cholQop,x-mu))
logdet = len(x) * np.log(nugget) - detQ(cholQ) + detQ(cholQop)

print("logdet_chol",round(logdet,4))
print("quad_chol",round(quad,4))
print("like_chol",round(-0.5 * (quad + logdet),4))

logdet_chol -1261.3613
quad_chol 1061.1439
like_chol 100.1087
In [18]:
a = spdeRes.computeQuad()
print("-> Relative difference quadratic =",round((a-quad)/quad,6))

-> Relative difference quadratic = 0.000255
In [19]:
pcm = spdeRes.getPrecisionKriging()
a = detQ(cholQop)
b = pcm.computeLogDetOp(1,0)
print("log_det_op_chol",round(a,4))
print("log_det_op_api",round(b,4))
print("-> Relative difference =",round((b-a)/a, 4))

log_det_op_chol 12296.9389
log_det_op_api 11965.4372
-> Relative difference = -0.027
In [20]:
a = detQ(cholQ)
b = pcm.computeLogDetQ(10,0)
print("log_det_Q_chol",round(a,4))
print("log_det_Q_api",round(b,4))
print("-> Relative difference =",round((b-a)/a,4))

log_det_Q_chol 11255.7151
log_det_Q_api 11255.7151
-> Relative difference = 0.0
In [21]:
a = -0.5 * (quad + logdet)
b = spdeRes.computeLogLike(100,0)
print("likelihood api",round(a,4))
print("likelihood_chol",round(b,4))
print("-> Relative Difference =",round(((b-a)/a),4))

likelihood api 100.1087
likelihood_chol 21.9765
-> Relative Difference = -0.7805