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()
spde.init(model,grid,None,gl.ESPDECalcMode.SIMUNONCOND)
spde.compute()
Grid query¶
In [6]:
iuid = spde.query(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.query(dat)
Kriging¶
In [9]:
spdeRes = gl.SPDE()
spdeRes.init(model,grid,dat,gl.ESPDECalcMode.KRIGING,mesh)
spdeRes.compute()
In [10]:
iuid = spdeRes.query(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 1070.1094 like_chol 95.6259
In [18]:
a = spdeRes.computeQuad()
print("-> Relative difference quadratic =",round((a-quad)/quad,6))
-> Relative difference quadratic = 0.00245
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 12198.9601 -> Relative difference = -0.008
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 95.6259 likelihood_chol 10.5533 -> Relative Difference = -0.8896