Stochastic Partial Derivative Equations¶

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,1)
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.getPrecisionOpCs().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)
print("Variance =",round(np.var(gridExt["simuManual"][np.where(gridExt["NewSel"]==1)]),4))

Variance = 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["*estim"],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"]
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_quad = spdeRes.computeQuad()
print("-> Relative difference quadratic =",round(100*(a_quad-quad)/quad,2),"%")

-> Relative difference quadratic = 0.03 %
In [19]:
pcm = spdeRes.getPrecisionKrig()
a_op = detQ(cholQop)
b_op = pcm.computeLogDetOp(1,0)
print("log_det_op_chol",round(a_op,4))
print("log_det_op_api",round(b_op,4))
print("-> Relative difference =",round(100*(b_op-a_op)/a_op, 2),"%")

log_det_op_chol 12296.9389
log_det_op_api 12296.9389
-> Relative difference = 0.0 %
In [20]:
a_one = detQ(cholQ)
b_one = pcm.computeLogDetQ(10,0)
print("log_det_Q_chol",round(a_one,4))
print("log_det_Q_api",round(b_one,4))
print("-> Relative difference =",round(100*(b_one-a_one)/a_one,2),"%")

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(100*(b-a)/a,2),"%")

likelihood api 100.1087
likelihood_chol 99.974
-> Relative Difference = -0.13 %
In [ ]: