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 gstlearn.document as gdoc
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

gdoc.setNoScroll()
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 [ ]: